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After  a  period  of  fast  development,  procedures  for  the  probabilistic 
quantification  of  seismic  hazard  have  reached  a  state  of  relative  maturity. 
Standard  methods  usually  differ  only  by  secondary  characteristics.  In  spite 
of  this  operational  stability,  estimates  of  seismic  hazard  that  are  produced  by 
different  algorithms  depend  on  and  are  often  very  sensitive  to  the  form  of  the 
mathematical  model  and  the  choice  of  some  of  the  parameters.  It  is  important 
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to  recognize  this  sensitivity  in  order  to  avoid  placing  excessive  confidence  in 
numerical  estimates  of  hazard,  especially  in  the  case  of  rare  events.  Errors 
may  originate  from  conscious  simplification  of  the  mathematical  representation 
(e.g.,  from  the  use  of  a  Poisson  model  of  earthquake  occurrences  as  opposed  to 
a  physically  more  attractive  point  process  with  spatial  and  temporal  memory), 
from  limited  statistical  information,  or  from  faulty  judgment  (expert  errors). 
The  purpose  of  this  report  is  to  review  different  types  and  sources  of  errors 
and  to  assess  their  impact  on  calculated  seismic  hazard.  In  the  available 
literature,  no  generally  accepted  procedure  exists  to  systematically  deal  with 
uncertainty  on  seismic  hazard  and,  most  frequently  in  practice,  elusive  argu-  , 
ments  of  conservatism  are  used  to  avoid  direct  confrontation  with  the  problem. 

Scientific  studies  are  sparse  and  fragmentary;  they  typically  focus  on 
narrow  aspects  of  the  problem,  use  different  statistical  techniques,  and  report 
results  in  different  formats.  Accordingly,  an  effort  is  made  here  to  coordinate 
these  results  and  provide  a  unified  picture.^ 

A  first  and  fundamental  problem,  which  ih  the  past  has  received  only 
minimal  attention,  is  to  give  a  definition  of  objective  or  "true"  hazard 
against  which  various  estimates  can  be  evaluated.  One  such  definition  is  pro¬ 
posed  here,  based  on  the  notion  of  hazard  as  relative  frequency  in  time.  Three 
classes  of  errors  are  then  introduced  and  separately  analyzed:  1.  Errors  due 
to  uncertainty  of  model  parameters,  such  as  the  geometrical  configuration  of 
the  earthquake  sources,  the  parameters  of  the  magnitude-frequency  relationship 
(including  the  upper  bound),  and  those  of  the  attenuation  function.  A 
recently  proposed  analytical  procedure  for  the  analysis  of  such  errors  is  re¬ 
viewed  and  results  from  sensitivity  studies  are  presented.  2.  Errors  due  to 
the  use  of  a  wrong  form  of  model.  These  errors  result  in  biases  (hence,  their 
statistical  magnitude  does  not  depend  on  the  amount  of  information  available) 
and  may  be  very  large,  especially  for  rare  events.  3.  Finally,  errors  of  in¬ 
terpretation  and  modeling  by  seismologist  experts.  Case  studies  are  reviewed 
and  mathematical  models  for  errors  of  this  last  type  are  proposed. 
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1 .  INTRODUCTION 

For  centuries,  engineering  design  has  been  guided  by  the  continuity  of 
successful  practice  and  professional  judgement.  A  first  element  of  objec¬ 
tive  rationality  has  been  provided  by  the  principles  of  mechanics,  the 
application  of  which  has  become  easier  through  the  use  of  fast  computing 
devices.  A  second  movement  toward  rationalization,  based  on  the  formal 
treatment  of  quantifying  uncertainty,  is  now  under  way.  By  itself,  none 
of  the  three  elements--mechanics,  computers,  and  quantitative  analysis  of 
uncertainty— can  produce  intelligent  engineering  decisions,  but  their 
combination  seems  to  be  both  necessary  and  sufficient.  Benefits  from 
probabilistic  treatment  of  uncertainty  are  especially  large  in  problem 
areas,  such  as  earthquake  engineering,  where  solutions  by  traditional 
means  require  a  great  deal  of  experienced  judgment  and  intuition. 

1 . 1  The  Problem  of  Accuracy  in  Seismic  Hazard  Analysis 

Wherever  a  mathematical  model,  deterministic  or  probabilistic,  is 
used  to  represent  a  complicated  physical  phenomenon,  the  problem  of  de¬ 
termining  its  accuracy  arises.  It  is  clear  what  accuracy  means  in  the 
case  of  deterministic  models  (one  should  just  compare  deterministic  pre¬ 
dictions  with  actual  observations);  this  is  not  so  for  probabilistic 
models.  For  example,  there  are  statisticians  who  consider  probability  in 
seismic  hazard  statements  as  quantified  personal  uncertainty  about  future 
earthquake  events.  From  this  premise  they  conclude  that  if  an  individual 
faithfully  quantifies  his  own  state  of  uncertainty,  then  calculated  hazard 
is  necessarily  exact,  at  least  for  that  individual.  This  argument  is 
built  around  the  so-called  subjective,  degree-of- belief  (d.o.b.)  interpre¬ 
tation  of  probability.  This  Interpretation  will  be  discussed  at  length  in 
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Sec,  2  and  conclusions  will  be  drawn  that  it  is  inappropriate  or  at  least 
insufficient  in  the  context  of  seismic  hazard.  One  should  notice,  however, 
that  the  issue  of  accuracy  is  legitimate  within  the  d.o.b.  interpretation 
for  two  reasons:  (1)  because  it  is  typical  to  use  simple  mathematical  ap¬ 
proximations  of  complicated  states  of  uncertainty  and  (2)  because  there  is 
ample  evidence  that  humans  make  mistakes  in  quantifying  their  own  state  of 
uncertainty. 

It  will  be  proposed  in  Sec.  2.1  that  probability  be  given  a  frequentist, 
time-average  interpretation  in  the  quantification  of  seismic  hazard;  e.g., 
that  the  probability  that  the  peak  ground  acceleration  at  a  given  site  dur¬ 
ing  the  next  50  years  will  exceed  0.2  g  be  interpreted  as  the  relative 
frequency  with  which  the  event  occurs  over  periods  of  50  years.  Notice 
that  in  this  interpretation,  the  next  time  interval  is  replaced  with  the 
generic  time  interval.  The  same  interpretation  has  two  major  advantages 
over  the  d.o.b.  interpretation:  (1)  it  is  objective  and  (2)  it  applies 
if  the  physical  earthquake  process  is  regarded  as  deterministic.  It  also 
makes  clear  what  is  meant  by  an  accurate  model:  an  accurate  model  is  one 
which  produces  probabilities  of  seismic  events  close  to  the  actual  relative 
frequencies  in  time. 

Accuracy  of  a  procedure  for  the  estimation  of  seismic  hazard  is  not 
by  itself  a  sufficient  measure  of  engineering  value.  ("Seismic  hazard"  is 
used  here  to  denote  a  suitable  probabilistic  characterization  of  future 
ground  motions  at  the  site,  often  simply  in  terms  of  the  mean  rate  at  which 
peak  motion  amplitudes  at  the  site  exceed  given  values.)  What  is  of  greater 
interest  for  engineering  purposes  is  "seismic  risk,"  i.e.,  the  probabilistic 
characterization  of  consequences  from  exposure  of  engineering  facilities  to 
seismic  hazard.  For  example,  it  may  be  of  interest  to  know  the  yearly  prob¬ 
ability  of  overtopping  of  a  dam  due  to  earthquake  events.  If  the  dam  is 


completely  unsafe  with  respect  to  very  severe  ground  motions,  then  in¬ 
accuracy  of  seismic  hazard  in  that  range  is  inconsequential. 


A  second  reason  why  accuracy  is  not  sufficient  is  that  simpler  pro¬ 
cedures  are  often  preferred  to  more  precise  but  complicated  methods. 
Decision-theoretic  approaches  to  engineering  model  selection  have  been 
proposed  (e.g.,  Grigoriu,  Veneziano,  and  Cornell,  1979)  and  they  account  for 
nrhe!  errors,  their  consequences,  and  the  cost  of  developing  and  using  the 
model.  It  would  be  interesting  to  examine  the  accuracy  of  seismic  hazard 
methods  from  this  wider  decision-theoretic  perspective,  but  studies  on  this 
subject  are  almost  nonexistent.  For  some  types  of  engineering  systems,  in¬ 
formation  is  available  on  the  dependence  of  risk  on  seismic  hazard  and  also 
on  the  sensitivity  of  risk  to  errors  on  seismic  hazard.  This  information  is 
reviewed  in  Sec.  3,  as  it  is  useful  to  evaluate  the  practical  relevance  of 
hazard  estimation  errors. 

1.2  Seismic  Hazard  Procedures  and  Sources  of  Errors 

Seismic  hazard  procedures  or  models,  M,  are  sequences  of  analytical 
or  numerical  operations  that  produce  estimates  of  earthquake  hazard  from 
suitable  seismological  information.  Denote  by  z  the  available  information 

A 

and  by  H  the  hazard  estimate  from  z  and  from  a  given  procedure  M.  The 
form  of  H  varies  depending  on  the  purpose  of  seismic  hazard  analysis:  for 
example,  H  might  consist  of  a  random  sequence  of  earthquake  times  and  in¬ 
tensities  at  the  site.  In  other  cases  H  is  simply  a  "hazard  function," 
e.g.,  a  function  A  (a)  which,  for  each  given  peak  ground  acceleration  a, 
gives  the  mean  rate  at  which  a  is  exceeded  at  the  site.  In  other  cases 

/S 

still,  H  may  characterize  seismic  threat  at  several  or  even  at  all  points 
of  a  geographical  region.  Estimates  of  this  last  type  are  needed  to 


! 
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characterize  the  vulnerability  of  spatially  extended  systems  (networks, 
lifelines,  groups  of  facilities),  to  establish  insurance  risk,  or  to  study 
the  socio-economic  impact  of  earthquakes. 

Although  procedures  of  different  types  have  been  proposed,  a  certain 
degree  of  standardization  has  been  achieved  in  recent  years  in  favor  of 
models,  the  elements  of  which  are:  a  set  of  "seismic  sources"  near  the  site 
of  interest,  a  random  process  model  of  time,  location,  and  size  of  earth¬ 
quakes  from  each  source,  and  an  attenuation  function  that  relates  ground 
motion  amplitude  at  the  site  to  earthquake  size  and  distance  from  the  energy 
source  (Cornell,  1968;  Der  Kiureghian  and  Ang,  1977).  In  the  present  study, 

A 

only  procedures  of  this  standard  type  are  considered  and  H  is  assumed  to  be 
a  hazard  function,  e.g.,  of  the  type  previously  denoted  by  A(a). 

In  order  to  understand  in  which  ways  H  may  be  in  error,  it  is  useful 
to  distinguish  among  the  following: 

My  =  true  (actual)  but  unknown  model .  This  is  a  model  which  produces 
probabilities  of  hazard- related  events  equal  to  the  associated 
relative  frequencies  in  time.  The  associated  true  hazard  is 
denoted  by  Hy. 

®  =  class  of  (simple)  seismic  hazard  models  to  which,  for  the  purpose 
of  obtaining  H,  the  true  model  My  is  assumed  to  belong.  For 
example,  in  =  {M(x,  ej}  =  class  of  Poisson  occurrence  models  with 
constant  mean  activity  rate  X  in  space  and  time  and  independent 
identically  distributed  (I.I.D.)  earthquake  sizes  Y-j ,  Y2,  ...  . 

The  Y.j  have  a  given  probability  distribution  function  Fy(y,  £) 
with  e  an  unknown  vector  of  parameters.  Hence,  in  this  case 
is  a  family  of  Poisson  models  parameterized  by  X  and  e_. 
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M  =  estimate  of  My  obtained  from  given  information,  under  the  assump¬ 
tion  that  My  belongs  to  Hi.  For  example,  M  =  M(X,  6j  in  which 

/>S 

X  and  e_  are  point  estimates  of  X  and  e_.  The  hazard  estimate  H 
corresponds  to  the  model  M. 

The  empirical  validation  of  H  or  the  evaluation  of  estimation  errors 
by  comparing  predictions  from  H  with  observations  from  Hy  is  impractical 
because  seismic  events  of  engineering  interest  are  too  rare.  More  realistic 
means  of  accuracy  evaluation  are  as  follows: 

If  indeed  My  belongs  to  in ,  then  one  can  study  properties  of  the  random 
estimation  error  H-Hy  by  means  of  Bayesian  analysis  or  by  methods  of  clas¬ 
sical  statistics.  In  the  former  case,  the  properties  of  the  error  depend 
on  the  distribution  of  My  in  HI;  in  the  second  case,  the  same  properties 
depend  on  which  of  the  elements  of  HI  is  the  true  model .  Of  course,  in  both 
cases  the  error  distribution  depends  on  the  type  and  amount  of  information 
available. 

A  simpler  but  less  precise  method  for  this  type  of  error  analysis 
consists  of  calculating  the  sensitivity  of  H  to  the  choice  ofTfl.;  for  ex¬ 
ample,  sensitivity  with  respect  to  the  form  of  the  attenuation  and  size- 
frequency  functions. 

A 

If,  contrary  to  the  assumption  in  the  calculation  of  H,  My  does  not 
belong  to  HI  (e.g.,  if  the  actual  sequence  of  times,  locations,  and  earth¬ 
quake  sizes  has  memory,  and  Poisson  processes  with  independent  sizes  are 
the  only  models  in  Hi),  then  severe  bias  in  hazard  estimation  may  result. 

If  bias  is  small  over  a  range  of  erroneous  modeling  assumptions,  then  the 
hazard  procedure  is  said  to  be  robust  over  that  range.  Robustness  is  an 
important  property  because  in  practice,  My  is  always  known  not  to  belong 
to  Hi. 
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Hazard  estimation  errors  when  My  is  an  uncertain  model  in  Hi  are  con¬ 
sidered  in  Sec.  4;  those  when  My  does  not  belong  to  Si  are  studied  in  Sec.  5 

Reliance  on  statistical  analysis  only  with  historical  data  treated  as 
a  random  sample  may  lead  to  exceedingly  large  uncertainty  on  some  critical 
parameters  (typically  source  geometry  and  upper  bound  earthquake  size)  and 
hence  on  calculated  hazard.  One  way  to  reduce  this  variability  is  to  use 
parameter  and  uncertainty  estimates  by  seismologists  based  on  statistical 
as  well  as  geophysical  and  geological  information.  Expert  opinion  is  es¬ 
pecially  valuable  in  regions  of  low  seismicity  for  which  historical  samples 
are  small.  Whenever  statistical  analysis  is  possible,  experts  should  be 
provided  with  results  from  it.  Thus,  use  of  expert  opinion  is  not  to  be 
viewed  as  an  alternative  to  statistical  analysis,  but  rather  as  a  comple¬ 
menting  technique. 

Expert  opinion  is  not  itself  infallible,  as  demonstrated  by  the  vari¬ 
ability  of  estimates  from  different  experts.  It  is  therefore  important  to 
understand  any  model  errors  from  using  this  approach  and  to  devise  methods 
to  reduce  these  errors  by  combining  estimates  from  different  experts. 

This  is  the  subject  of  Sec.  6. 

Unfortunately,  errors  in  seismic  hazard  analyses  have  not  been  studied 
systematically;  contributions  are  sparse  and  fragmentary,  each  focussing 
on  a  narrow  aspect  of  the  problem.  A  variety  of  techniques  is  used  to  ex¬ 
press  results  in  different  forms.  In  spite  of  the  effort  to  give  unity  to 
the  subject,  voids  and  inconsistencies  still  remain.  There  is  no  accepted 
procedure  to  rationally  and  systematically  deal  with  uncertainty  on  seismic 
hazard  and,  frequently  in  practice,  elusive  arguments  of  "conservatism"  are 
used  to  avoid  direct  confrontation  with  the  problem. 
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I  Additional  work  is  needed  on  each  topic  in  this  report  and  a  special 

I  effort  should  be  made  to  unify  the  language  and  coordinate  the  results. 

I  Practical  guidelines  on  how  to  deal  with  non-frequentist  uncertainty  would 

I  be  of  much  value.  One  possibility,  which  is  advocated  by  many,  is  to 

regulate  safety  in  terms  of  events,  the  probability  of  which  is  high  enough 
to  be  accurately  estimated  by  available  procedures.  For  example,  one  might 

-3 

base  seismic  design  on  the  earthquake  intensity  with  a  return  period  of  10 
years,  as  opposed  to  a  return  period  of  106  years.  Of  course,  higher  pro¬ 
tection  against  damage  should  be  provided  under  the  less  intense  motion. 

It  is  hoped  that  considerations  in  this  study,  especially  those  on 
the  interpretation  of  seismic  hazard,  the  classification  of  estimation 
errors,  and  the  treatment  of  expert  opinion,  wi 1 1  prove  useful  to  future 


endeavors. 
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2.  FREQUENTIST  SEISMIC  HAZARD  AND  UNCERTAINTY  DUE  TO  ERRORS 

The  argument  that  uncertainty  expressed  by  seismic  hazard  is  the 
product  of  ignorance  and  that  seismic  hazard  probability  therefore  should 
be  regarded  as  a  measure  of  strength  of  belief  is  unsatisfactory  in  many 
respects  because  degrees  of  belief  (1)  do  not  necessarily  satisfy  the 
axioms  of  probability  theory,  (2)  are  vaguely  defined  and  variable  from 
individual  to  individual,  and  (3)  are  difficult  to  quantify.  Due  to  (1), 
the  mathematical  theory  of  probability  does  not  strictly  apply.  Diffi¬ 
culty  (2)  makes  it  questionable  whether  important  engineering  decisions 
should  be  based  on  essentially  subjective  hazard  assessments,  and  (3)  in¬ 
dicates  that  implementation  faces  problems  of  inference. 

In  Sec.  2.1,  rephrasing  of  seismic  hazard  statements  is  proposed, 
thus  making  it  possible  to  interpret  probability  as  relative  frequency 
in  time.  Such  an  interpretation  overcomes  the  above  difficulties  and,  con¬ 
trary  to  the  ensemble-frequentist  definition,  does  not  refer  to  a  statisti¬ 
cal  population,  which  in  the  present  case  obviously  does  not  exist. 

Not  all  uncertainties  that  contribute  to  seismic  hazard  are  of  the 
same  type  or  admit  the  time-average  frequentist  interpretation;  for  example, 
uncertainty  about  seismicity  parameters  does  not.  In  this  case,  how  should 
one  combine  probabilities  with  different  meanings  and  what  interpretation 
should  one  give  to  the  results?  We  shall  address  these  questions  in  Sec.  2.2 
by  analyzing  three  different  ways  of  dealing  with  multiple  interpretations: 
(1)  avoid  combination,  (2)  combine  probabilities  and  accept  a  weaker  interp¬ 
retation  of  the  results,  and  (3)  avoid  explicit  treatment  of  non-frequentist 
uncertainty  by  using  so-called  predictive  analysis.  The  second  alternative 
relies  on  the  hierarchical  order  that  exists  among  some  probability  inter¬ 
pretations.  Much  of  the  contents  of  this  section  is  from  current  research 
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by  Veneziano  and  Chung  (1980). 

2.1  The  Time-Average  (t.a.)  Interpretation  of  Probability 

In  the  case  of  seismic  hazard,  the  interpretation  of  probability  as 
degree  of  belief  can  be  avoided.  Consider,  for  example,  the  sequence  of 
times  at  which  strong  earthquakes  occur.  It  is  typically  found  that  these 
earthquakes  take  place  without  significant  regularity  or  clustering  in 
time,  in  a  way  that  has  been  effectively  described  by  Poisson  point  pro¬ 
cesses  (the  specific  form  of  the  process  has  no  relevance  in  the  present 
argument).  Although  the  Poisson  model--as  any  random  point  process--cannot 
describe  the  deterministic  physical  sequence  exactly,  there  are  at  least  two 
ways  in  which  it  may  be  correct: 

1.  It  may  be  correct  in  the  sense  of  degree  of  belief  (d.o.b. 

correct)  if  it  exactly  quantifies  uncertainty  about  (not 
intrinsic  in)  actual  seismicity  or,  more  important  here, 

2.  It  may  be  correct  in  the  frequentist  sense  that  probabilities 

calculated  from  the  Poisson  model  coincide  numerically  with 
limits  of  relative  frequency  in  time.  In  this  case,  the 
Poisson  model  is  t.a.  correct.  As  the  definition  makes 
obvious,  t.a.  models  have  the  form  of  stationary  random 
processes. 

For  the  second  interpretation  to  be  possible,  the  physical  sequence  of 
earthquakes  must  possess  some  degree  of  regularity,  so  that  limits  of 
relative  frequency  in  time  exist.  This  condition  does  not  imply  time- 
invariance  of  the  process  that  generates  the  natural  sequence;  geophysical 
laws,  earthquake  mechanism,  tectonic  stress  or  strain  rates  might  vary  in 
time  and  cause  clustering  of  the  earthquake  events,  and  yet  a  stationary 
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t.a.  process  may  exist.  Of  course,  in  this  case  the  process  must  have 
memory  and  not  be  one  of  the  Poisson  type.  Limits  of  relative  frequency 
in  time  satisfy  the  axioms  of  probability  theory. 

The  t.a.  interpretation  of  probability  is  not  new;  it  was  one  of  the 
first  interpretations  proposed  in  the  physical  sciences  (Hofstetter,  1964). 

For  example,  it  was  used  to  represent  the  erratic  variation  of  voltage  at 
the  terminals  of  a  resistor  and  of  velocity  at  a  given  point  in  turbulent 
air  flow.  Its  spatial  analogue  (probability  as  limit  of  relative  frequency 
in  space)  was  exploited  by  Matheron  (1971)  to  justify  certain  geostatistical 
procedures  of  mineral  exploration  and  estimation. 

Let  x(t)  denote  the  physical  function  of  interest,  e.g.,  x(t)  =  average 
shear  stress  on  the  plane  of  a  fault  at  time  t  .  It  may  not  be  easy  to 
prove  that  for  a  given  function  x(t)  an  associated  t.a.  process  X(t)  exists, 
but  this  is  of  little  consequence  here.  In  application,  x(t)  is  not  known 
and  assumptions  are  made  directly  on  the  random  process  X(t). 

In  this  context,  one  faces  the  opposite  problem,  i.e.,  to  establish 
whether  a  given  random  function  X ( t )  is  a  legitimate  t.a.  process  (whether 
or  not  a  function  x(t)  exists  that  has  X(t)  as  its  associated  t.a.  process). 

It  is  clear  that  t.a.  processes  do  exist.  For  example,  all  processes 
X(t)  that  are  ergodic  in  the  sense  that  mean  squares  are  legitimate 
t.a.  processes. 

It  is  also  clear  that  there  are  infinitely  many  deterministic  functions 
that  are  compatible  with  any  legitimate  t.a.  process  X(t);  let  x(t)  be  one 
such  function.  Then,  translates  of  x(t)  and  functions  that  differ  from  x(t) 
only  inside  a  finite  interval  of  time  still  have  X(t)  as  their  t.a.  process. 
More  importantly,  and  as  shown  by  the  theory  of  ergodic  processes,  the  family 
of  deterministic  functions  that  share  the  same  t.a.  process  includes  functions 


of  most  diverse  form.  Nor  should  one  only  think  of  deterministic  functions 
that  are  realizations  of  stationary  processes.  Because  of  time-averaging, 
functions  with  nonstationary  character  also  belong  to  this  family.  This 
indeterminacy  of  x(t),  given  X(t),  is  the  mechanism  through  which  t.a. 
probability  represents  uncertainty. 

The  point  should  be  stressed  that,  although  there  are  mathematical 
links  between  the  theory  of  ergodic  processes  and  the  partial  character¬ 
ization  of  an  unknown  function  through  its  t.a.  process,  there  are  also  im¬ 
portant  differences.  In  the  first  case,  one  views  the  physical  function 
itself  as  a  random  process,  whereas  in  the  second  case,  one  deals  with  the 
t.a.  properties  of  a  class  of  deterministic  functions. 

The  t.a.  interpretation  requires  some  modification  of  seismic  hazard 
statements.  Suppose  that  seismic  hazard  is  expressed  in  terms  of  the 
probability  distribution  of  =  peak  ground  acceleration  at  a  given  site 
during  the  next  t  years.  The  ensemble-frequentist  interpretation  is  not 
applicable  in  this  case  because  a  statistical  population  of  seismic  his¬ 
tories  at  the  site  cannot  be  defined.  In  t.a.  analysis,  the  definition  of 
A.  must  be  changed  and  fa. (a)  =  P[At  <  a]  interpreted  as  the  relative  fre- 
quency  in  time  with  which  a  is  exceeded  during  intervals  of  t  years. 

Although  specific  reference  to  the  next  time  interval  is  lost  due  to  time¬ 
averaging,  one  may  argue  that  this  is  a  very  reasonable  interpretation  of 
seismic  hazard,  if  degrees  of  belief  are  to  be  avoided.  Reference  to 
selected  intervals  of  time  can  still  be  made  by  way  of  conditional  t.a. 
analysis,  as  shown  next. 

Making  seismic  hazard  predictions  on  the  basis  of  the  stationary  process 
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X(t)  alone  may  not  be  appropriate  if  one  has  specific  information  on  the 
present  state  of  the  physical  system  (e.g.,  from  knowledge  of  the  recent 
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seismic  history  at  the  site).  Information  of  this  type  can  be  incorporated 
into  t.a.  hazard  statements  by  way  of  conditional  probabilities.  With  A 
and  B  as  any  two  events,  the  conditional  t.a.  probability  P[A|B]  is  defined  as 
the  limit  of  relative  frequency  in  time  with  which  event  A  occurs,  given 
that  event  B  occurs.  It  satisfies  P[A|B]  =  P[AnB]/P[B],  in  which  all  proba¬ 
bilities  are  t.a.  probabilities.  In  application  to  seismic  hazard,  B  rep¬ 
resents  information  on  the  physical  function  x(t)  other  than  that  already 
contained  in  the  random  process  X(t).  It  often  happens  that  event  B  occurs 
with  zero  frequency  in  time,  so  that  P[B]  =  P[ATIB]  =  0,  hence  the  previous 
equation  for  P[A]B]  cannot  be  used.  In  this  case  one  must  replace  proba¬ 
bilities  with  probability  densities  that  have  the  t.a.  meaning  of  densities 
of  relative  frequency  in  time. 

In  the  estimation  of  seismic  hazard,  all  available  information  on  the 
physical  function  x ( t )  should  be  used,  not  only  to  infer  properties  of  the 
associated  t.a.  process  X ( t ) ,  but  also  to  update  seismic  hazard  by  calcu¬ 
lating  conditional  t.a.  probabilities.  This  is  the  mechanism  through  which 
t.a.  analysis  can  express  time  dependence.  It  is  only  when  no  information 
on  x(t)  is  available--other  than  knowledge  of  the  process  X(t)--that  seismic 
hazard  predictions  are  independent  of  current  time. 

The  t.a.  interpretation  makes  it  possible  to  give  frequentist  meaning 
to  probabilistic  statements  of  seismic  hazard  and  suggests  itself  as  a 
measure  of  rational  d.o.b.,  whenever  applicable.  This  is  an  extension  of 
the  notion  that  whenever  the  ensemble  relative-frequency  interpretation  is 
possible,  the  beliefs  of  a  rational  individual  should  be  consistent  with 
frequentist  probabilities. 
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2.2  Combination  of  Frequentist  and  Non-Frequentist  Uncertainty 

The  frequentist  interpretation  of  seismic  hazard  is  conditional  on  a 
given  t.a.  process  representation  of  seismicity,  symbolized  here  by  the 
random  process  X(t).  However,  X(t)  is  not  always  known.  For  example,  one 
may  know  that  X(t)  is  a  homogeneous  Poisson  impulse  process  with  random 
event  sizes  Yi  at  the  times  t^  of  earthquake  occurrences. 


X(t)  =  £  Yi  6(t-ti) 

but  parameters  of  the  process  (the  mean  rate  of  events  and  the  marginal 
distribution  of  the  I.I.D.  variables  Y^)  may  not  be  exactly  known.  In  other 
cases,  the  form  of  the  process  X(t)  may  also  be  unknown. 

Uncertainty  about  the  form  and  the  parameters  of  the  t.a.  representa¬ 
tion  of  seismicity  cannot  be  measured  by  probabilities  with  frequentist  in¬ 
terpretation.  This  fact  causes  some  difficulty;  if  probability  is  to  be 
used,  then  it  becomes  necessary  to  allow  for  subjective  interpretations  and 
the  problem  arises  whether  or  not  different  interpretations  "mix,"  i.e., 
whether  or  not  in  the  course  of  analysis,  d.o.b.  probabilities  can  be  mean¬ 
ingfully  combined  with  frequentist  probabilities.  How  should  results  be 
interpreted?  This  problem,  which  has  important  implications  on  the  pro¬ 
cedure  of  seismic  hazard  analysis  and  on  the  format  and  numerical  value  of 
hazard  results,  is  currently  the  source  of  much  controversy. 

Advocates  of  the  d.o.b.  interpretation  observe  that  measures  of  strength 
of  belief  are  numerically  compatible  with  frequentist  probabilities  whenever 
a  frequentist  interpretation  is  possible.  Therefore,  all  probabilities-- 
including  those  that  result  from  mixing  operations--can  be  interpreted  in 
the  same  d.o.b.  sense. 
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Other  statisticians  disagree,  either  because  they  object  to  quanti¬ 
fying  the  d.o.b.  through  probability,  or  because  they  cannot  attribute 
any  meaning  to  mixtures.  At  most,  they  accept  statements  of  the  type 


'The  d.o.b. 

that  A(a)  £  A  is  pdob(x)" 

(a) 

The  d.o.b. 

that  jpfr[At  £  a]  <  p] 

r  <*  pdob<P>" 

(b) 

In  the  first  equation,  A(a)  is  the  mean  rate  at  which  peak  ground  accelera¬ 
tion  at  the  site  exceeds  the  value  a.  Subscripts  fr  and  d.o.b.  denote  fre¬ 
quent  ist  (for  example,  t.a.)  and  d.o.b.  interpretations,  respectively. 

In  case  (b),  •|pfr[At  £  a]  £  pj  is  viewed  as  an  event  to  which  a  d.o.b. 
probability  is  assigned.  Hence  in  this  case— that  of  multiple  interpreta¬ 
tions  without  mixing--it  is  possible  to  talk  of  the  "probability  of  a  proba 
bility,"  in  the  sense  of  the  d.o.b.  in  the  relative  frequency"  of  an  event. 

With  reference  to  Eq.  2.1b,  practical  interest  is  in  P[Afc  £  a].  If  one 
decides  that  mixing  is  possible,  e.g.,  because  frequentist  probabil ities 
also  have  the  meaning  of  d.o.b.,  then  P[At  <  a]  can  be  found  from  the  Total 
Probability  Theorem,  which  gives 


P[At5 


p  dPdob{p) 


(2.2) 


We  have  therefore  identified  three  different  attitudes  toward  the 
use  and  mixing  of  probability  interpretations: 
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1.  Different  interpretations  are  allowed  and  probabilities  can 
be  mixed  (Eq.  2.2). 

2.  Different  interpretations  may  coexist  without  mixing  (Eq.  2.1). 

3.  Only  frequentist  interpretations  are  allowed. 

In  order  to  focus  on  conceptual  issues,  we  have  been  purposely  schematic 
in  classifying  probability  interpretations  a*.  L-ier  "frequentist"  or 
d.o.b."  Other  interpretations  exist  afM  ,s(v;etimes  subtle)  variants 
have  been  proposed  with  the  frequentist  *  .fe.  frameworks  (Fine,  1973; 
Lucas,  1970).  In  the  case  of  more  interpretations,  generalization 

of  Eq.  2.1  (especially,  the  second  statement)  produces  very  confusing 
language,  and  option  2  should  be  avoided. 

Option  1  (mixing)  seems  to  be  available  only  for  probabilities  with 
identical  interpretations,  and  since  a  number  of  interpretations  may  be 
acceptable  for  the  same  probability,  one  should  first  determine  whether 
it  is  possible  to  have  a  single  common  interpretation.  A  hierarchical 
order  exists  for  some  probabilities;  for  example,  frequentist  interpreta¬ 
tions  are  stronger  than  d.o.b.  interpretations,  since  for  rational  in¬ 
dividuals  the  former  imply  the  latter,  whereas  the  reverse  is  not  true. 
Similarly,  the  rational  d.o.b.  interpretation  is  stronger  than  the  personal 
d.o.b.  interpretation.  Ordering  is  not  possible  in  all  cases,  e.g.,  there 
is  no  ordering  between  the  ensemble  and  the  t.a.  frequentist  interpretations 

It  may  be  that  common  interpretations  are  too  weak  to  be  useful.  For 
example,  one  might  interpret  all  probabilities  as  personal  d.o.b. ,  thus 
losing  objectivity  and  having  to  face  problems  mentioned  earlier  in  this 
section.  Before  mixing  probabilities,  one  should  therefore  look  for  the 
strongest  possible  common  interpretation.  Also,  in  the  context  of  seismic 
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hazard,  many  seismologists  and  engineers  seem  to  prefer  statements  of  the 
type  illustrated  by  Eq.  2.1  to  the  alternative  of  mixing  probabilities  and 
producing  results  with  weaker  d.o.b.  interpretation. 

Since  mixing  (Eq.  2.2)  reduces  the  strength  of  interpretation  and 
multiplicity  of  meanings  (Eq.  2.1)  produces  unappealing  statements,  it  would 
be  valuable  to  develop  procedures  that  combine  uncertainties  of  different 
types  and  still  retain  a  strong  interpretation  of  the  results  (Option  3).  This 
would  seem  to  be  an  impossible  task  in  the  light  of  what  has  been  just 
said.  However,  some  such  procedures  are  available  which,  for  specific 
problems,  account  for  parameter  uncertainty  using  frequentist  analysis  ex¬ 
clusively.  Results  are  in  the  form  of  prediction  regions  of  given  proba¬ 
bility  content  and  are  obtained  by  the  method  exemplified  next. 

Consider  for  simplicity  a  stationary,  independent  sequence  of  variables 
Xi  (i=l,2,...),  the  marginal  CDF  of  which,  Fx(x),  depends  on  an  unknown 
vector  of  parameters  £.  This  might  be  the  sequence  of  earthquake  intensi¬ 
ties  at  the  site.  Statistical  information  on  £  is  obtained  by  observing 
n  variables  of  the  sequence,  and  we  want  to  characterize  uncertainty  on 
the  next  unobserved  variable  Xp+^ .  One  can  do  this  by  calculating  "pre¬ 
diction"  sets  on  the  real  line  with  given  probability  content.  A  predic¬ 
tion  set  of  P-content  is  a  random  set  D  =  D ( X^ ,  ...,  Xp)  such  that  for  any 
given  £,  D  has  probability  content  P,  i.e.. 


•  Xr 


P[D(X1,...,Xn)|0]  = 


all  x 


.*r 


n 

n 

i=l 


dF 


xle 


(x,-) 


.  f  dFx|i(x)  =  P 
D(x1,.'.,xn) 


(2.3) 
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Notice  that  Eq.  2.3  treats  the  quantities  Xj,...,X  as  random  variables 
and  that  uncertainty  in  _0  is  not  being  quantified.  In  particular,  one 
may  constrain  D  to  be  a  one-sided  interval  of  the  type  (-°°,  5(Xp...,X  ,P)]. 
Then  £(X-j ,. . .  ,Xn,P)  is  the  random  P-fractile  of  the  so-called  predictive 
distribution  of  Xn+^ .  As  the  sample  size  n  increases,  uncertainty  in  6_  de¬ 
creases  and  prediction  regions  tend  to  become  smaller.  It  should  be  em¬ 
phasized  that  frequentist  prediction  sets  may  not  exist  (the  problem  is 
analogous  to  that  of  finding  confidence  intervals  for  unknown  parameters). 

The  literature  is  rich  in  results  on  frequentist  prediction  regions 
for  univariate  and  multivariate  sequences  and  for  simple  and  multiple 
(simultaneous)  prediction  (Proshan,  1953;  Aitchison,  1964;  Aitchison  and 
Sculthorpe,  1965;  Aitchison  and  Dunsmore,  1975;  Chew,  1968;  Guttman, 

1970).  For  a  review,  see  Veneziano  (1974).  In  the  case  of  multiple  pre¬ 
diction  (Chew,  1968),  D(X.| , . . .  ,Xn)  is  a  random  set  in  the  sample  space  of 
m  unobserved  variables  from  the  sequence,  Xn+1 ,  ....  Xn+m. 

Frequentist  prediction  regions  are  often  numerically  identical  or  very 
close  to  Bayesian  prediction  regions  if,  for  the  latter,  0  has  (essentially) 
non informative  prior  distribution. 

The  definition  of  the  frequentist  prediction  region  in  Eq.  2.3  applies 
also  to  probabilities  with  t.a.  frequentist  meaning. 

2.3  Representation  of  Seismic  Hazard 

It  is  proposed  that  an  objective  probabilistic  measure  of  seismic  hazard 
be  defined  as  the  frequency  in  time  of  relevant  earthquake  events.  The  pur¬ 
pose  of  hazard  analysis  is  to  estimate  this  objective  probability  from  given 
information. 

If  information  is  limited,  then  estimation  is  imperfect  and,  whenever 
possible,  one  should  obtain  "predictive"  frequentist  statements  of  hazard 


A 


20 


which  incorporate  uncertainty  on  seismicity  parameters  without  recourse 
to  d.o.b.  interpretations.  If  frequentist  predictive  statements  are 
impossible,  then  an  alternative  is  to  quantify  estimation  uncertainty 
through  d.o.b.  probability  on  true  hazard.  That  is,  hazard  probabilities 
are  themselves  considered  as  uncertain  quantities  with  non-frequentist  dis¬ 
tributions.  Bayesian  analysis  produces  probabilities  that  correspond  to 
the  mean  value  of  these  distributions.  The  mean  value  is  not  necessarily 
the  best  estimate  of  seismic  hazard;  the  very  fact  that  many  engineers  do 
not  make  decisions  only  on  the  basis  of  Bayesian  hazard  indicates  that  other 
characteristics  of  the  non-frequentist  distribution  of  hazard  are  also  im¬ 
portant.  When  viewed  from  a  "no-mixing"  perspective,  the  Bayesian  approach 
is  correct  but  incomplete;  it  is  still  useful  to  quantify  the  bias  of  pro¬ 
cedures  that  neglect  estimation  errors.  In  this  sense  we  shall  use  Bayesian 
analysis  in  Sec.  3.2.2.  Knowledge  of  the  non-frequentist  distribution  of 
hazard  is  important  to  assess  the  impact  of  future  research  and  data-gathering 
efforts  aimed  at  improving  knowledge  of  seismicity  and  procedures  of  hazard 
analysis. 
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3.  RISK  CONSEQUENCES  OF  ERRORS 

Typical  seismic  hazard  procedures  generate  single  (point)  estimates 
of  frequentist  hazard,  rossibly  wrong  assumptions  (e.g.,  spatial  and 
temporal  independence  of  earthquake  events)  and  uncertainty  on  the  value 
of  parameters  (e.g.,  of  the  magnitude-frequency  and  intensity-attenuation 
functions)  are  sources  of  estimation  errors.  Before  quantifying  these 
errors  in  Secs.  4,  5,  and  6,  we  shall  look  at  the  influence  that  uncer¬ 
tainty  on  seismic  hazard  has  on  engineering  risk.  Should  one  find  that 
risk  is  insensitive  to  this  uncertainty,  the  importance  of  errors  in 
hazard  analysis  would  be  much  reduced. 

The  present  task  is  complicated  by  two  facts:  (1)  non-frequentist 
uncertainty  on  hazard  depends  on  the  seismic  province  and  site  in  question, 
and  (2)  seismic  risk  depends  on  seismic  resistance,  which  is  itself  uncer¬ 
tain  and  variable  from  facility  to  facility.  Nevertheless,  some  informa¬ 
tion  about  typical  uncertainty  on  hazard  and  variability  of  seismic  resistance 
is  available,  thus  allowing  exploratory  treatment  to  the  problem.  This  can  be 
done  by: 

1.  Fitting  parametric  seismic  hazard  models  and  resistance  distributions 
to  available  data  (Sec.  3.1). 

2.  Calculating  seismic  risk  with  parameters  either  fixed  to  their  mean 
values  (thus  ignoring  estimation  errors,  Sec.  3.2.1)  or  treated  as 
random  variables  (thus  including  the  effect  of  uncertain  estimation 
errors.  Sec.  3.2.2). 

3. 1  Simple  Models  of  Seismic  Hazard  and  Seismic  Resistance 


For  the  purpose  of  calculating  engineering  risk  (e.g.,  the  probability 
or  the  mean  rate  of  "failure"  events),  seismic  hazard  and  resistance  should 


22 


be  described  in  terms  of  the  same  ground  motion  intensity  parameter,  Y. 

For  example, 

Y  =  in  a,  where  a  =  peak  ground  acceleration,  or 

Y  =  in  v,  where  v  =  peak  ground  velocity,  or 

(3.1) 

Y  =  Jin  Q(u),  £),  where  Q(o),  g)  =  maximum  pseudo-velocity  of  an 

oscillator  with  natural  frequency  to  and  damping  ratio  £,  or 

Y  =  I,  where  I  =  Modified  Mercalli  (MM)  intensity 

Models  of  Seismic  Hazard.  Much  experience  has  been  accumulated  on 
hazard  functions  in  terms  of  the  above  intensity  measures.  Some  such  func¬ 
tions  obtained  for  Boston  under  different  modeling  assumptions  are  shown  in 
Fig.  3.1  as  relationships  between  MM  intensity  I  and  the  mean  rate  A(I)  at 
which  events  of  larger  local  size  occur  (Cornell  and  Merz,  1975). 

When  plotted  on  semilog  paper  ( Jin  A  versus  Y),  hazard  functions  are 
typically  well  approximated  by  straight  lines  in  the  low-to-medium  intensity 
range  but  decay  faster  at  higher  intensities.  Let R denote  the  resistance  of  the 

statistical  population  of  facilities  under  consideration,  expressed  in  the 

? 

same  units  as  Y,  and  denote  by  mR  and  oR  the  mean  value  and  Vc-.riance  of  R, 
respectively.  Then  failure  occurs  if  Y  >  R  or  equivalently  i!  Y1  >  R',  in 
terms  of  the  normalized  variables 


R'  =  (R-mR)/aR 
Y'  =  (Y-mR)/oR 


(3.2) 


A  few  simple  analytical  hazard  functions  of  the  type  A’(y’)  =  mean  rate  at 
which  normalized  site  intensity  y'  is  exceeded  are  given  in  Table  3.1  and 
plotted  in  Fig.  3.2  for  selected  parameter  values.  The  associated  hazard 


(FIRM  GROUND)  M  M  INTENSITY  IN  BOSTON 


FIGURE  3.1  Influence  of  different  assumptions  on 
seismic  hazard  in  Boston  (from  Cornell 
and  Merz,  1975). 
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CASE 

- mm 

MODEL 

PARAMETERS 

a 

Linear 

>•  e-6'*' 

b 

Truncated  linear 

V  e-6'y',  y'  $  yi’ 
o,  y  >  y,' 

A'.g'.y^ 

c 

Truncated  linear 
dA ' (y ' )/dy ’ 

x'(e*0  y  -e"e  yl  \y’<y1 1 
0,  y'^' 

A’ ,B’ ,y1 • 

d 

Truncated  linear 
d2A' (y1 )/dy*2 

X •  [(e-B,y ’ -e_8  yl ' )  ( 

-  8'(y-|  '-y')e'6  yl  J.y'^y, ' 

0,  y^y/ 

A*.  B'OV 

e 

Quadratic 

x.  g-a'y^-B'y’ 

X* ,a* ,8’ 

TABLE  3.1 


Simple  normalized  hazard  models. 


FIGURE  3.2  Untruncated  and  tru 
models  In  Table  3.1 
yj  is  assumed  zero. 
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functions  in  terms  of  abnormalized  intensity,  A(y),  can  be  found  from 

A(y)  =  A'((y-mR)/oR)  (3.3) 

For  example,  in  the  case  of  the  first  model  in  Table  3.1,  Eq.  3.3  gives 

A(y)  =  A  e'By  (3.4) 

e'mR 

in  which  A  =  A'e  and  6  =  B'/aR. 

The  hazard  functions  of  Table  3.1  are  flexible  enough  to  accurately 

fit  nearly  all  practical  cases.  For  example,  the  hazard  curves  in  Fig.  3.1 

are  well  approximated  by  "linear"  models  with  slope  parameter  B  *  1.1  (curves 

UB  12  and  RANDOM  12)  or  6  =  1.77  (curve  CA  12),  or  by  "truncated  linear" 

models  with  6  between  1.7  and  2.0  (remaining  curves).  More  sophisticated 

models  in  Table  3.1  would  provide  even  better  fits.  They  would  also  be 

conceptually  more  satisfactory  because,  as  observed  by  many  authors, 

there  should  be  a  physical  upper  limit  to  intensity;  however, 

the  truncated  linear  model  unrealistically  associates  a  nonzero  mean  rate 
-BV 

(of  magnitude  A'e  ')  with  this  limit. 

Models  of  Seismic  Resistance.  Errors  in  the  hazard  function  may  have 
substantially  different  consequences  on  risk,  depending  on  the  resistance 
characteristics  of  the  exposed  facility.  Safety  engineers  have  been  es¬ 
pecially  concerned  with  estimation  accuracy  at  high  ’evels  of  intensity, 
hence  for  return  periods  1/A  that  require  extrapolation  beyond  the  duration 
of  historical  records.  This  concern  is  justified  by  the  increase  of  un¬ 
certainty  on  unA  given  y  (on  y  given  a)  with  increasing  y  (with  decreasing  A) . 
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However,  rare  high  intensity  events  may  not  be  of  determining  impor¬ 
tance  in  the  calculation  of  risk. 

Consider  a  population  of  engineering  facilities  (e.g.,  rockfill  dams 
or  high-rise  buildings)  the  seismic  resistance  R  of  which  has  probability 
density  function  fR{y).  The  associated  cumulative  distribution  function, 

ry 

FR(y)  =  fR(r)dr,  is  sometimes  referred  to  as  the  "fraqil ity  curve. "  It 

—co 

gives  the  probability  of  failure  (the  fraction  of  population  that  fails)  if 
an  earthquake  occurs  with  site  intensity  y.  Also,  denote  by  X (y )  =  dA(y)/dy 
the  mean  rate  density  of  events  with  site  intensity  exactly  y,  and  let  seis¬ 
mic  risk  be  measured  by  the  mean  failure  rate  X^..  For  example,  a  value 
Xf  =  10"2  failure/year  means  that  if  facilities  from  the  population  were 
exposed  to  probabilistically  independent  but  otherwise  identical  seismic 
environments,  the  average  life  of  the  generic  facility  would  be  E[T]  = 

1/A.f  =  100  years.  Under  the  assumption  of  Poisson  events,  the  failure 

-Xft 

probability  of  the  generic  facility  in  t  years  would  be  P^(t)  =  1-e 

The  definition  of  R  as  the  resistance  of  the  generic  facility  of  a  sta¬ 
tistical  population  makes  risk  results  meaningful  in  a  frequentist  sense.  If, 
nowever,  only  one  facility  exists,  then  fR  expresses  non-frequentist  uncer¬ 
tainty  on  the  resistance  of  that  facility  and  calculated  risk  can  have 

_2 

only  d.o.b.  interpretation.  For  example,  X^  =  10  failures/year  must  be 
taken  to  mean  that  the  d.o.b.  distribution  of  time  to  failure  has  a  mean 
value  of  100  years.  In  order  to  avoid  statements  of  this  type,  which  are 
often  obscure,  uncertainty  on  seismic  resistance  will  always  be  given  fre¬ 
quentist  interpretation  here. 

The  mean  failure  rate  can  be  calculated  by  any  of  the  following 


expressions: 
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fR(y)A(y)dy 

(a) 

FR(y)*(y)dy 

(b) 

fR.ty'jA'ly'Jdy' 

(c) 

F^y’Jx'fy'Jdy1 

(d) 

(3.5) 


in  which  R1  is  the  standardized  resistance  in  Eq.  3.2  and  X1  is  the  deriva¬ 
tive  of  A' . 

Errors  may  be  present  not  only  in  the  estimated  hazard  (functions  A 
and  X)  but  also  in  the  estimated  resistance  distribution  (functions  F  and 
f).  According  to  Eqs.  3.5a  and  3.5b,  X^.  is  a  weighted  average  of  FR,  fR, 

A,  or  X  with  weight  functions  X,  A,  fR,  and  FR,  respectively.  Hence,  for 
example,  in  order  to  evaluate  the  seismic  risk  consequences  of  errors  in 
A(y)  for  different  intensity  values  y,  one  should  look  at  the  corresponding 
weight  function,  fR:  where  the  latter  function  is  maximum  (at  the  mode  of 
R)  errors  in  A  are  most  critical.  This  is  actually  true  for  errors  of 
fixed  magnitude,  but  in  the  case  of  A,  errors  tend  to  be  larger  for  smaller 
intensities  (notice  that  interest  here  is  in  A,  not  in  tn  A).  One  con¬ 
cludes,  then,  that  the  most  critical  region  for  accuracy  of  A  is  below  the 
mode  of  R. 

It  is  also  of  interest  to  look  at  the  whole  integrands  of  Eqs.  3.5a 
and  3.5b,  as  they  vary  with  y.  The  integrand  of  the  first  equation, 
fR(y)My).  gives  the  contribution  to  the  failure  rate  from  facilities  with 
different  resistance,  whereas  the  integrand  of  the  second  equation. 


Fp(y)X(y),  gives  the  contribution  to  from  earthquakes  of  different 
intensity.  The  first  function  is  maximum  below  the  mode  of  R;  this  is 
usually  true  also  for  the  second  function. 

The  early  belief  that  uncertainty  on  seismic  resistance  contributes 
little  to  overall  risk  and  that  even  models  with  deterministic  resistance 
produce  valuable  results  has  been  seriously  questioned.  Analysis  in  Sec. 

3.2  also  indicates  that  variability  of  R  may  be  critical,  especially  when 
site  intensity  has  an  upper  bound  or  statistical  uncertainty  on  resistance 
distribution  parameters  is  included.  It  appears  that  for  some  engineering 
facilities  a  substantial  fraction  of  risk  comes  from  medium-size  earthquakes 
that,  although  associated  with  small  failure  probabilities,  are  much  more 
frequent  than  large  and  more  destructive  events. 

At  the  present  time,  the  probabilistic  prediction  of  seismic  resis¬ 
tance  by  explicit  analytical  means  (e.g.,  random  vibration)  is  unreliable 
for  all  but  the  simplest  systems.  Repeated  deterministic  analyses, using 
simulated  ground  motions  conditional  on  site  intensity  and  simulated  sys¬ 
tem  behavior  and  resistance  parameters,  are  more  general  but  also  very 
expensive  procedures.  For  some  types  of  facilities  such  as  earth  dams, 
building  frames,  and  major  structural  components  of  nuclear  power  plants, 
simplified  analytical  procedures  (Banon  and  Veneziano,  1981;  Hasselman  and 
Simonian,  1980)  that  may  prove  feasible  and  sufficiently  accurate  are  now 
under  development.  In  many  cases,  however,  the  best  source  of  information 
is  historical  seismic  performance.  Much  of  this  information  has  been  col¬ 
lected  in  terms  of  the  damage  ratio  (DR),  which  is  the  cost  of  repair 
divided  by  the  cost  of  replacement.  Data  on  wooden  frame  buildings  and  on 
ordinary  or  reinforced  masonry  construction  has  been  available  for  some 
time.  In  the  last  few  years,  damage  statistics  on  high-rise  buildings, 
large  dams,  and  important  engineering  facilities  have  also  been  gathered, 
in  part  (e.g.,  in  the  nuclear  industry)  through  large-scale  testing  (see 
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Proceedings  of  SMiRT  Conferences:  Berlin,  1971,  1973;  London,  1975; 

San  Francisco,  1977).  For  high-rise  buildings,  reliable  information 
has  been  collected  from  the  1971  San  Fernando  earthquake.  The  most 
extensive  survey  (Whitman,  1973)  documents  368  buildings  of  five  stories 
or  more,  classified  by  age,  structural  material,  and  height.  Fig.  3.3 
displays  mean  damage  ratio  (MDR)  data  on  high-rise  buildings  from  the  San 
Fernando  as  well  as  other  earthquakes.  The  only  classification  here  is 
by  UBC  zone  where  all  heights  and  construction  materials  (steel  and 
concrete)  are  lumped  together.  Based  in  part  on  statistical  data,  in 
part  on  professional  judgement,  curves  relating  MDR  to  MM  intensity  have 
been  proposed  by  several  authors.  Those  in  Fig.  3.4  for  high-rise  build¬ 
ings  are  from  Whitman  (1973).  Other  estimates  can  be  found  in  Mann  (1974), 
Whitman  and  Hong  (1973),  and  Benjamin  (1974).  For  most  cases  and  over  a 
wide  range  of  intensities,  the  expected  log  damage  ratio  varies  almost 
linearly  with  I  and  can  be  represented  as 

E[£n  DR | I ]  =  ad  +  bpl  (3.6) 

with  constants  ap  and  bp  that  depend  on  the  type  of  facility. 

Some  information  is  also  available  on  the  dispersion  of  the  condi¬ 
tional  damage  ratio,  DR | I .  Benjamin  (1974)  found  that  for  broad  classes 
of  buildings,  the  damage  data  for  given  intensity  are  fitted  well  by  either 
lognormal  or  gamma  distributions.  If  the  lognormal  model  is  used,  then 
(fcn  DR | I )  has  normal  distribution.  The  same  author  found  that  the  vari¬ 
ance  of  (£n  DR | I )  is  approximately  constant  with  I.  Damage  statistics  in 
Whitman  (1973)  confirm  that  °£noR | i  is  not  sensitive  to  I  and  is  compati¬ 
ble  with  a  normal  distribution  of  (£n  DR | I ) . 
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FIGURE  3.4  Mean  damage  ratio  for  high-rise  buildings 
designed  under  different  UBC  zone  require¬ 
ments.  S  denotes  "superzone"  with  twice 
the  lateral  force  for  Zone  3  (from  Whitman, 
1973). 
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In  later  analye;s,  it  is  not  the  value  of  0£nQR|j 
rather  the  ratio 


that  counts,  but 


6 


D 


£n  DR  1 1 


(3.7) 


Table  3.2  collects  some  statistics  and  subjective  evaluations  of  the 
parameters  bp  and  8p.  The  values  from  Newmark  (1974)  and  Vanmarcke  (1971) 
are  not  strictly  compatible  with  those  from  Benjamin  (1974)  and  Whitman 
(1973)  because  the  former  refer  to  given  peak  ground  acceleration,  a,  in¬ 
stead  of  given  MM  intensity.  For  example,  Newmark  suggests  values  of 
a&n( response)/9  for  orc|inary  buildings  and  nuclear  reactor  structures  and 
equipment  (parameter  BETA  in  his  Table  3).  If  the  level  of  response  is 
proportional  to  a  and  a  varies  by  a  factor  of  2  per  unit  of  I,  then  the 
response  varies  also  by  a  factor  of  2  per  unit  of  I  and  estimates  of  Bp  in 
Eq.  3.7  can  be  found  from 


6 


0 


in  2 

°sin(response)/a 


(3.8) 


This  last  relationship  has  been  used  to  calculate  values  of  Bp  in  Table  3.2 
from  estimates  of  a*n(response) (a. 

In  making  a  model  of  probabilistic  resistance,  the  following  assump¬ 
tions  are  introduced: 

1.  "Failure"  occurs  when  the  damage  ratio  reaches  a  critical  value  D*; 

2.  The  conditional  expected  £n  damage  ratio  is  a  linear  function  of  I, 
of  the  type  in  Eq.  3.6;  and 

3.  The  conditional  distribution  of  sin  given  I  is  normal,  with  con- 

2 

stant  variance  Op. 


Benjamin  (1974) 

upper  bound 
lower  bound 

0.484 

0.347 

mm 

Whitman  (1973) 

Post-1947  Buildings 

San  Fernando,  1=6 

=1.13 

=1.88 

San  Fernando,  1=7 

*1.13 

=2.13 

San  Fernando,  1=7.5 

=0.91 

=1.90 

Newmark  (1974) 

Nuclear  Reactor  Structure 

■I 9 

Nuclear  Reactor  Equipment 

n§ 

Vanmarcke  (1971) 

1=6.7 

=  0.68 

1.24 

1=7.8 

=0.59 

1.53 

TABLE  3.2  Values  of  bQ  and  in 
Eqs.  3.6  and  3.7. 
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It  follows  from  1,  2,  and  3  that  the  probability  of  failure  from  an  earth¬ 
quake  of  intensity  I  is 


Pf(I)  =  0  [(D*-aD-bDI)/oD] 


(3.8) 


in  which  <f>  is  the  standard  normal  CDF.  From  the  fact  that  Pf(I)  is  also 

the  probability  that  R  <  I,  one  concludes  that,  when  expressed  in  units  of 

2 

MM  intensity,  R  has  normal  distribution  N(mR,  oR),  with  parameters 


mR  "  "aD^bD 

2  _  2 /h2  _  0-2 
°R  “  °D'bD  6D 


(3.9) 


Also,  in  terms  of  normalized  intensity  Y'  =  I*  =  (I-mR)/oR,  resistance  R' 
has  standard  normal  distribution. 

Results  on  the  distribution  type  for  R  are  not  expected  to  change  if 
ground  motion  intensity  at  the  site  is  measured  in  terms  of  any  other  param¬ 
eter  Y  in  Eq.  3.1 . 

3.2  Calculation  of  Seismic  Risk 

Given  the  models  of  seismic  hazard  in  Table  3.1  and  given  normality 
of  the  distribution  of  R,  the  mean  failure  rate  Xf  is  conveniently  calcu¬ 
lated  from  Eq.  3.5c.  Results  should  first  be  obtained  when  all  uncertainty 
is  of  the  frequentist  type  (when  the  parameters  of  A'  and  fR  are  known) 
and  then  consider  a  few  cases  in  which  parameter  uncertainty  is  included. 

3.2.1  Cases  with  Known  Seismic  Hazard 

Combinations  of  R'  ~  N(0,  1)  with  any  of  the  functions  A'  in  Table  3.1 
are  considered  first.  Sensitivity  to  the  distribution  of  R'  is  then  checked 
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by  assuming  that  resistance  has  gamma  distribution  and  that  A'  has  untrun¬ 
cated  linear  form.  In  all  these  cases,  A^  is  found  by  closed-form  integra¬ 
tion  of  Eq.  3.5c. 
a.  Linear  Hazard 

,  co  -6  'y  * -y 1  2/2  B,2/2 

Af  =  -.•=  e  dy'  =  A'e  (3.10) 

f  ^2"  J.co 


For  hazard  expressed  in  units  of  MM  intensity  and  using  Table  3.2  for  Bq, 
one  finds  that  the  parameter  6'  =  6/Bq  (B  is  the  slope  of  the  seismic  hazard 
curve)  varies  between  0.6  and  1.2,  with  a  typical  value  of  0.9.  The  constant 

A'  in  Eq.  3.10  is  the  mean  rate  at  which  earthquakes  occur  with  intensity  in 

d'2/o 

excess  of  the  mean  resistance;  hence,  the  quantity  y  =  ep  can  be  viewed 
as  a  penalty  factor  for  uncertain  resistance;  it  increases  with  oR  (=  1/6q) 
and  with  6  and  is  1  if  either  parameter  is  zero.  In  the  present  model, 
typical  values  of  y  are  in  the  range  1.2  to  2.1,  but  much  higher  values  may 
result  from  using  nonlinear  A'  models  or  from  including  statistical  uncer¬ 
tainty  on  Op  into  the  analysis, 
b.  Truncated  Linear  Hazard 

In  this  case  A^  is  given  by 


Xf  = 


yl 

~f 

/2rr  '  -OO 


dy. 


=  A' 


B'2/2 


4>(y1  '+6') 


(3.11) 


in  which 
S'  2/2 

e  i^y-j  +  6 1 ) ,  which  includes  the  effect  of  truncation,  is  plotted  in 


yj  =  (y-|"mR)/0R  is  normalized  upper  bound  intensity.  The  factor 


37 


Fig.  3.5  as  a  function  of  yj  for  selected  values  of  6'.  For  yj  >1, 
the  mean  failure  rate  in  Eq.  3.11  is  essentially  the  same  as  for  the  un¬ 
truncated  case.  In  order  to  reduce  Xf  to  1/2  the  value  for  untruncated 
hazard,  it  must  obviously  be  that  y-|  =  m^  -  o^B,  i.e.,  the  upper  bound 
intensity  at  the  site  must  be  about  one  standard  deviation  of  resistance 
below  the  mean  resistance. 

An  example  of  application  to  Boston  is  given  in  Fig.  3.6  and  Table  3.3. 
The  figure  shows  five  truncated  or  untruncated  linear  approximations  to  the 
hazard  curves  of  Fig.  3.1.  For  the  case  of  resistance  R  (in  units  of  MM 

o 

intensity)  with  normal  distribution  N(8,  (0.8)  ),  the  parameters  8'  =  Bop, 
X',  and  y^  =  (y-]-mp)/Op  for  the  five  approximations  are  given  in  Table  3.3, 
columns  2,  3,  and  4,  respectively.  Column  5  gives  the  associated  failure 
rate.  Consideration  of  the  modeling  assumptions  behind  each  case  indicates 
that  truncation  of  the  site  intensity-frequency  law  has  relatively  small 
effect  on  X ^  unless  truncation  is  significantly  below  the  mean  resistance 
(compare  cases  2  and  3).  However,  truncation  of  the  epicentral  intensity- 
frequency  relationship  is  more  important  because  it  also  reduces  the  value 
of  X1.  Compare  case  1  (truncation  at  I  =  XII  for  all  the  sources)  with  cases 
3  and  5  in  which  upper  bound  intensities  between  6.3  and  7.3  were  assumed 
for  seismic  sources  near  Boston.  The  increase  of  B'  due  to  truncation  of 
epicentral  intensity  has  only  a  minor  balancing  effect. 

Integration  of  Eq.  3.5c  for  cases  c  and  d  in  Table  3.1  gives  the  fol¬ 
lowing  failure  rates. 

c.  Truncated  Linear  First  Derivative  of  A 


Xf  = 


X' 


eB'2/2  $(y1,+B')  - 


-B1 


*1 


* (y, ' ) 


(3.12) 


MMI  at  Boston  (Firm  Ground) 


FIGURE  3.6  Linear  and  truncated  linear  approximations 
to  the  hazard  functions  in  Fig.  3.1. 


HAZARD 
CURVE 
(Fig.  3.6) 

8'  =  Bor 

1 

0.88 

5.9*  10' 


8.7*  10 


8.6*  10' 


2.9  *  10 


8.6  *  10' 


2.7  *  10' 


2.6*  10 


»  9.1*10 


1.58  2.6  *  1  O'5  -0.625  j  5.3*  10'5 


TABLE  3.3  Parameters  and  failure  rate  for  the 
five  hazard  curves  in  Fig.  3.6. 
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d.  Truncated  Linear  Second  Derivative  of  A 


=  x'  |e3'2/2  4>(y-j  '+&')  «  e  P  yi  [(l+e'y-, '  ‘ )  +  3,'4>(y1  ’ )]]■ 


-e’y 


(3.13) 


in  which  <p  is  the  standard  normal  probability  density  function.  The 
parameter  X1  in  the  last  two  equations  is  the  mean  rate  of  failure  for 
R  =  mR,  in  the  case  of  no  truncation,  i.e.,  y-|-*-°°. 

One  should  be  cautioned  against  directly  comparing  the  mean  crossing 
rates  of  the  last  three  models  because  when  using  them  to  approximate  actual 
hazard  functions,  one  would  conceivably  select  different  values  of  y^  and 
possibly  also  of  X'  and  e'. 

e.  Quadratic  Hazard 


X 


f 


21  f  dy 


X' 


(3.14) 


Drawbacks  of  the  quadratic  hazard  model  are  that  it  has  no  upper  bound  and 
that  it  violates  the  condition  that  A  should  be  a  nonincreasing  function  of 
y.  Results  using  this  model  are  typically  very  close  to  those  from  model  c. 

f.  Linear  Hazard  and  Gamma  Distribution  of  Resistance 

In  order  to  test  sensitivity  of  risk  to  the  shape  of  the  distribution 
of  R,  suppose  that  R'  has  shifted  gamma  distribution  with  zero  mean  and 
unit  variance,  hence  with  probability  density  function 
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0  , 


fR,(y')  = 


P[D(-y'+D)] 

r(D2) 


D  -1 


-D(y'+D) 


y1  <  -D 

y 1  >  -D 


(3.15) 


in  which  D  is  a  positive  constant.  Plots  of  this  function  are  shown  in 
Fig.  3.7  for  selected  values  of  D.  Notice  that  the  above  gamma  distribution 
reduces  to  a  shifted  exponential  for  D  =  1  and  to  the  standard  normal  for 
D-*°°.  For  A'  that  corresponds  to  the  untruncated  linear  model,  the  mean 
failure  rate  is 


f  -  I  "  e-rs  d,' 

r(D^)  J-D 


-  A  1  Pe'°  D 

"  A  e  U+B1" 


O' 


(3.16) 


The  ratio  between  this  rate  and  that  for  normal  resistance  distribu¬ 
tion  with  the  same  mean  value  and  variance  (see  Eq.  3.10)  is  given  by 


¥9”)  f  D  1°2  6*D-B'2/2 
^(normal )  B  [d+B '  J  e 


(3.17) 


This  ratio  is  shown  in  Fig.  3.8  as  a  function  of  B'  for  D  =  1,  2,  3,  5,  8. 
The  failure  rate  ratio  is  generally  smaller  than  1,  primarily  because  the 
gamma  density  vanishes  for  y'  <  -D.  This  feet  shows  the  importance  of  left 
truncation  of  the  resistance  distribution.  However,  if  one  excludes  the 
artificious  cases  when  D  <  1  and  for  typical  values  of  6'  (e.g..  S'  <  1.5), 
the  mean  failure  rate  is  insensitive  to  the  present  change  of  resistance 


FIGURE  3.8  Plots  of  the  mean  failure  rate 
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distribution.  Higher  sensitivity  is  found  with  respect  to  distributions 
that  allocate  more  probability  to  values  of  R  below  the  mean. 

3.2.2  Cases  with  Uncertain  Seismic  Hazard 

Non-fre quentist  uncertainty  on  the  mean  failure  rate  \f  is  contributed 
by: 

1.  Possible  errors  in  the  seismic  hazard  procedure  (M^  might  not  belong  to 
I;  see  Sec.  1.2).  One  might  call  these  "model  errors"; 

2.  Uncertainty  on  which  model  in  Hi  corresponds  to  the  true  model  fL,  given 
that  e  HI .  If  Hi  is  a  parametric  family  of  models,  then  this  source 
of  error  might  be  called  "parameter  uncertainty"; 

3.  Errors  in  the  form  of  the  distribution  of  R  and  on  the  resistance 

2 

parameters,  mR  and  oR. 

Except  for  very  few  cases,  the  distribution  of  induced  by  these 
possible  errors  and  uncertainties  must  be  found  through  straightforward 
but  tedious  numerical  integration,  which  requires  that  very  ma.iy  risk  cal¬ 
culations  be  made,  each  conditional  on  different  modeling  assumptions  or 
parameter  values.  The  fact  that  realistic  approximations  to  the  hazard 
function  lead  to  explicit  expressions  for  X^  suggests  an  alternative  fwo- 
step  procedure  which,  although  approximate,  is  much  more  expeditious.  First, 
one  quantifies  uncertainty  on  the  hazard  function  A (y )  due  to  possible 
errors  in  the  assumptions  and  to  uncertainty  on  the  parameters  of  the  seis¬ 
mic  hazard  procedure.  (Results  on  uncertainty  on  A  will  be  given  later  in 
Secs.  4,  5,  and  6.)  Then  one  relates  uncertainty  on  A(y)  to  the  probability 
distribution  or  at  least  to  the  expected  value  of  XR.  The  second  step  be¬ 
comes  particularly  simple  if  A(y)  has  the  form  of  model  a  in  Table  3.1,  with 
uncertain  parameters  X'  and  e'.  Uncertainty  on  the  parameters  of  the  resis¬ 
tance  distribution  can  also  be  considered.  The  results  that  follow  on  the  mean 
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of  Af  are  limited  to  this  case.  In  Bayesian  analysis,  one  is  interested 
only  in  the  probability  E[P^(t)]  and  for  small  risk  this  is  closely  approx¬ 
imated  by  E[A^]  •  t.  In  classical  analysis,  one  would  typically  be  interested 
in  the  full  distribution  of  A^  {or  of  P^(t)),  but  one  can  still  use  the  mean 

value  E[Af]  to  calculate  the  bias  of  estimates  that  do  not  account  for  un¬ 
certainty  on  seismic  hazard  or  on  resistance  distribution  parameters.  In 
some  cases  the  full  distribution  of  A^  is  obvious  from  its  dependence  on 
the  uncertain  parameters;  in  other  cases  it  can  be  found  in  approximation 
or  through  numerical  convolution. 

Under  the  assumption  that  A1  has  untruncated  exponential  form  (model 

2 

,  oR),  non-frequentist 
? 

uncertainty  is  limited  to  the  four  parameters  A',  B' ,  mR,  and  oR.  Results  on 
E[A^]  for  all  combinations  of  known/unknown  parameters  with  appropriate  dis¬ 
tributions  are  given  in  Veneziano  (1975);  those  for  known  resistance  distri¬ 
bution  are  reproduced  here  in  terms  of  the  ratio  y  between  E[A^]  and  the 
mean  failure  rate  when  the  uncertain  parameters  are  fixed  to  mean  values, 

Eq.  3.10.  Hence,  y  can  be  interpreted  as  a  bias  factor  due  to  parameter 

2 

uncertainty.  Let  mR  and  oR  be  known.  Then  for  uncertain  A1  or  B1  one  finds: 

1.  If  S'  is  known  and  A1  is  uncertain  with  lognormal  distribution  such  that 
in  A1  -  N  (m£nAlf  o\ny  )  ,  then 


a  in  Table  3.1)  and  that  R  has  normal  distribution  N(mR 


y  yj 

>.f  |  X 1  =exp(m. 


•  ex',(lW2! 


(3.18) 


2 

In  the  range  of  exp(cr£n^,)  from  0.4  to  1.1,  which  corresponds  to  a 
1-sigma  uncertainty  factor  on  A’  between  1.5  and  3,  the  ratio  y^,  var¬ 
ies  from  1.1  to  1.8.  In  this  case,  A^  has  itself  lognormal  distribu¬ 


tion. 
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n 


o 

2.  If  X'  is  known  and  8‘  has  normal  distribution  N(mDl,  cr;;,),  then  the 

P  P 

mean  failure  rate  ratio  Yg.  is 

V  =  171  =  n-o^)_1/2  exp[m2,  o2,/2(l-a2,)]  (3.19) 

For  realistic  distribution  parameters  (mD1  between  0.8  and  1.6  and  o„, 

8  P 

between  0.1  and  0.3),  y^,  varies  in  the  range  1.01  to  1.20,  with  the 
upper  limit  very  sensitive  to  the  assumed  maximum  value  of  oQl. 

P 

3.  If  X1  and  6'  are  independent  with  the  marginal  distributions  given 

above,  then  the  factor  is  simply  the  product  of  the  marginal 

factors,  y^,  =  Y^.y^i.  Most  often,  however,  the  assumption  of 

independence  is  inappropriate.  It  is  found  that  if  «.n  X'  and  6'  have 
correlation  coefficient  p,  then  the  factor  y^,  gl  can  be  decomposed 
into  three  partial  factors: 

V,e'  =  V{p)  •  V  ‘YP  (3-20) 

where  y^,  is  given  by  Eq.  3.19, 

Yx,(p)  =  exp[(l-p2)  <^nX./2] 

and 

Yp  =  exPtpOjmx^P^nx.+^g.OgJ^d-a2,)] 

For  p  =  0,  Eq.  3.20  reproduces  the  result  for  independent  parameters. 

If  p  =  1,  then  uncertainty  on  X'  is  totally  explained  by  6'  and 
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Ya.(1)  =  1;  if  in  addition  o£nX>  =  0,  then  X'  is  known  and  y^>  g.  =  Yg>- 

2 

Explicit  results  when  m^  and  are  uncertain  (with  conjugate  distributions 
in  random  sampling)  are  not  quoted  here.  Whereas  uncertainty  on  has 
moderate  consequences  on  E[X^],  uncertainty  on  may  be  very  important  be¬ 
cause  it  considerably  increases  the  probability  of  small  resistance  values  to 
the  point  that  for  exponential  hazard,  many  failures  may  result  from  ground 
motions  of  small  intensity.  For  detailed  results,  see  Veneziano  (1975). 

Although  limited  by  the  simple  forms  of  the  hazard  and  resistance  dis¬ 
tribution  functions,  results  of  this  section  should  be  kept  in  mind  when 
evaluating  the  consequences  of  errors  in  hazard  analysis;  in  particular, 
errors  that  affect  the  upper  bound  intensity  and  errors  in  the  seismicity 
rates  and  in  the  slope  of  intensity-frequency  relationships. 
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4.  ERRORS  FROM  UNCERTAINTY  ON  MODEL  PARAMETERS 

Error  analysis  in  engineering  seismic  risk  encompasses  two  prob¬ 
lems:  (1)  evaluation  of  non-frequentist  uncertainty  on  seismic  hazard, 
and  (2)  calculation  of  the  consequences  of  this  uncertainty  on  risk. 
Problems  of  the  second  type  have  been  examined  in  the  last  section,  T 
whereas  the  remainder  of  the  report  focuses  on  the  first  task.  This 
is  a  complicated  one  because  of  the  multiplicity  of  error  sources  and 
seismic  hazard  procedures.  It  is  subdivided  into  three  parts  and  each 
part  treated  in  a  separate  section. 

The  present  section  looks  at  the  variation  of  hazard  due  to  modeled 
uncertainty  on  or  to  reasonable  variations  of  the  parameters,  while 
Sec.  5  quantifies  the  consequences  of  unmodeled  errors.  Hence,  it  is 
assumed  here  that  the  true  model  My  is  of  the  type  hypothesized  in 
hazard  analysis  (that  My  e  1R )  whereas  the  opposite  assumption  is  made 
in  Sec.  5.  In  application,  there  is  a  tendency  to  rely  on  expert 
opinion  in  order  to  include  nonstatistical  information  and  reduce 
estimation  errors.  Expert-based  estimates  are  themselves  not  error- 
free.  Conceptual  and  mathematical  models  and  some  results  on  the 
variability  of  expert  opinion  will  be  presented  in  Sec.  6. 

Returning  to  the  objective  of  the  present  section,  the  aim  is  to 
quantify  the  uncertainty  on  hazard  (e.g.,  on  the  mean  exceedance  rate 
A(a))  due  to  the  uncertainty  of  the  form  and  parameters  of  the  seismicity 
model,  e.g.,  to  uncertainty  on  shape,  size,  and  location  of  earthquake 
sources,  on  the  activity  of  each  source  (magnitude  distribution  and 
upper  bound,  earthquake  occurrences  in  time  and  space),  on  the  form  and 
parameters  of  the  attenuation  function,  and  on  site  amplification  effects. 
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Use  of  probabilistic  analysis  to  propagate  uncertainty  through  hazard 
calculations  is  made  difficult  by  two  factors:  (1)  the  relationship 
between  the  uncertain  quantities  and  hazard  may  not  be  explicit  but 
rather  given  in  algorithmic  form  and  (2)  it  is  difficult  to  quantify 
uncertainty  about  some  aspects  of  the  model,  e.g.,  the  magnitude  dis¬ 
tribution  or  the  spatial  dependence  of  earthquake  epicenters.  As  a 
consequence,  one  must  resort  to  different  procedures  for  different 
parameters,  ranging  from  exact  analytical  treatment  to  numerical  sensi¬ 
tivity.  Powerful,  although  approximate,  analytical  techniques  have  been 
recently  developed  by  McGi'ire  and  Shedlock  (1980)  and  will  be  reviewed 
first  in  Sec.  4.1.  Results  from  parameter  sensitivity  analyses  given 
the  form  of  the  model  will  follow  in  Sec.  4.2.  Finally,  Sec.  4.3  will 
comment  on  the  effect  of  changing  basic  assumptions  in  the  hazard  model. 

4. 1  Analytical  Results 

Cases  of  parameter  uncertainty  that  lend  themselves  to  analytical 
treatment  become  fewer  as  complexity  of  the  model  increases.  The 
following  model,  simple  but  realistic,  has  been  studied  by  McGuire  and 
Shedlock  (1980). 

Earthquake  sources  are  independent.  They  generate  events  at  Poisson 
times,  with  truncated  exponential  distribution  of  magnitude.  Earthquakes 
originate  from  fault  ruptures  and,  given  magnitude,  rupture  length  has 
truncated  lognormal  distribution  (truncation  is  due  to  the  finite  length 
of  the  faults).  Peak  ground  motion  characteristics  at  the  site  are 
related  by  an  attenuation  fu».  •. ion  to  magnitude  and  shortest  distance  to 
the  rupture  segment.  Conditional  on  magnitude  and  distance,  they  too 
have  lognormal  distribution.  Uncertainty  on  hazard  at  the  site  is  studied 


by  treating  model  parar*  ^rs  either  as  known  constants,  or  as  mutually 
independent  random  variables. 

The  main  approximation  of  the  method  consists  of  treating  earth¬ 
quake  sources  as  points  in  space,  thus  replacing  the  random  distance 
to  the  site  with  a  given  constant  for  each  source.  Similarly,  single 
values  are  used  to  replace  magnitude  distributions. 

Seismic  hazard  is  given  in  terms  of  the  peak  ground  acceleration 

ap  that  corresponds  to  a  specified  annual  probability  of  exceedance, 

0 

V 

It  was  found  by  McGuire  and  Shedlock  that  the  distribution  of  aD 
due  to  parameter  uncertainty  has  a  mean  value  (=  Bayesian  estimate) 
practically  identical  to  the  value  calculated  by  fixing  all  parameters 
to  their  mean  value.  This  unbiased  property  of  conventional  hazard 
procedures  was  attributed  to  almost  linear  dependence  of  ap  on  the 
uncertain  parameters.  (But  notice  that  bias  may  be  large  if  there  are 
errors  in  model  assumptions ,  as  shown  in  Sec.  5,  or  if  hazard  is  nonlinear 
in  the  uncertain  parameters;  see  e.g.  Fig.  4.4). 

Sources  and  magnitude  of  non-frequentist  uncertainty  vary  from  case 
to  case.  In  an  application  to  the  San  Francisco  Bay  area,  McGuire  and 
Shedlock  considered  the  following  parameters  to  be  uncertain:  the  depth 
of  energy  release,  the  activity  rate  and  Richter  b  value  for  each  fault, 
and  the  mean  attenuation  relationship  from  magnitude  to  acceleration.  An 
additional  source  of  uncertainty  for  sites  in  the  central  Mississippi 
Valley  was  taken  to  be  the  exact  geographical  location  of  faults.  For 
each  active  fault,  the  method  that  these  authors  propose  to  calculate  the 
distribution  of  a0  due  to  parameter  uncertainty  is  as  follows. 
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First,  one  uses  standard  methods  of  hazard  analysis  to  calculate  ap 
as  a  function  of  P,  with  all  uncertain  parameters  and  relationships 
fixed  to  mean  values.  In  particular,  the  value  of  ap  for  P  =  Pq  is 
considered  to  be  the  mean  value  of  ap  when  parameter  uncertainty  is 
included. 

One  then  calculates  an  "effective  magnitude"  and  an  "effective 

distance"  Rg,  defined  as  the  mean  magnitude  and  mean  distance  of  events 

that  cause  peak  site  acceleration  larger  than  ap  ,  i.e., 

0 

M  =  I  M  *  P[A  >  ap  i M]  •  P[M]/P 
M  0 


Rp  =  l  R  *  P[A  >  ap  | R]  •  P[R]/P0 
e  R  K0  U 


(4.1) 


The  method  proceeds  by  calculating  sensitivity  factors  on  the  peak  site 

acceleration  ap  due  to  variations  of  the  uncertain  parameters.  Parameters 
0 

are  varied  one  at  a  time,  and  in  each  case,  one  utilizes  the  quantities  in 

Ea.  4.1  to  replace  lengthy  seismic  hazard  convolutions  with  closed-form 

calculations  using  point  sources  and  single  magnitude  values. 

For  each  parameter  combination  (parameters  are  assumed  to  have  discrete 

distribution)  one  finds  ap  by  using  the  best  estimate  of  ap  and  pre- 

K0  K0 
viously  calculated  sensitivity  factors.  One  also  finds  the  associated 

probability  as  the  product  of  individual  parameter  probabilities. 

The  final  result  is  a  discrete  distribution  of  ap.  In  the  case  of  a 

single  active  fault,  this  approximate  distribution  was  compared  with  the 

exact  distribution  from  repeated  seismic  hazard  analysis  for  each  parameter 

combination  and  found  to  be  very  accurate. 


I 
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Numerical  results  using  realistic  representations  of  uncertainty  are 

as  follows.  For  the  San  Francisco  Bay  area,  the  coefficient  of  variation 

of  an  V,  ,  ranges  from  less  than  0.2  at  sites  located  about 

a0.002 

50  km  from  the  main  faults  to  more  than  0.4  near  these  faults.  Inter¬ 
mediate  coefficients  of  variation  are  obtained  at  sites  far  away  from 
the  earthquake  sources.  The  reason  for  this  peculiar  dependence  of 

V  on  distance  R  from  the  high  intensity  region  is  that  a  large 
“0.002 

contribution  to  uncertainty  on  aD  comes  from  uncertainty  on  the  mean 

K0 

attenuation  function.  Data  to  which  attenuation  relationships  are 

fitted  are  abundant  at  distances  of  about  50  km;  and  for  these  distances, 

alternative  attenuation  functions  are  nearly  the  same.  For  sites  with 

R  «  50  km,  more  important  sources  of  hazard  uncertainty  ire  the  activity 

rate  and  the  Richter  b  value,  whereas  close  to  the  faults,  another 

important  source  of  variation  of  aD  is  the  depth  of  energy  release. 

K0 

In  the  calculation  of  V  for  the  central  Mississippi  Valley, 

“0.002 

uncertainty  on  the  actual  location  of  the  faults  was  also  modeled.  In 
this  case,  location  of  the  earthquake  sources  and  depth  of  energy  release 
were  found  to  be  the  main  sources  of  hazard  uncertainty,  especially  at 
sites  close  to  the  faults.  Away  from  the  epicentral  region,  uncertainty 
on  attenuation  becomes  relatively  more  important.  The  largest  values  of 

V  in  excess  of  0.4  were  found  in  the  regions  of  higher  activity. 
“0.002 

4.2  Sensitivity  Analyses 

Quite  likely,  methods  of  error  analysis  of  the  semianalytical  type 
described  in  the  previous  section  or  of  the  type  by  enumeration  will 
become  of  common  use  in  the  near  future.  They  allow  one  to  represent 
hazard  not  in  terms  of  a  single  frequency- intensity  curve,  but  as  a 
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continuous  or  discrete  family  of  curves,  each  associated  with  an  exceedance 
probability.  This  representation  of  seismic  hazard  is  being  advocated  by 
several  authors,  e.g.  Kaplan  (1980)  and  Kaplan  and  Garrick  (1980),  as  a 
way  to  express  confidence  in  the  results  of  standard  procedures  of  hazard 
analysis.  The  same  representation  corresponds  to  the  no-mixing  option  of 
Sec.  2.2;  it  has  the  conceptual  appeal  of  avoiding  the  combination  of 
probabilities  with  different  interpretations  and  of  separating  objective 
hazard  from  subjective  or  non-frequentist  uncertainty.  Its  main  draw¬ 
backs  are  (1)  the  nonavailability  of  practical  instruments  of  computation 
(the  algorithm  by  McGuire  and  Shedlock  is  the  first  of  this  kind  to  be  pro¬ 
posed)  and  (2)  the  added  complication  and  expense  of  analysis.  Also,  as 
observed  by  McGuire  and  Shedlock,  the  representation  of  seismic  hazard  is 
useful  when  making  decisions  about  research  and  data  gathering  efforts, 
but  it  is  not  clear  how  one  should  utilize  it  when  dealing  with  the  safety 
of  given  facilities.  In  the  latter  context,  it  seems  that  only  the  average 
hazard  counts.  If  this  is  so  and  if,  as  found  by  McGuire  and  Shedlock,  the 
average  hazard  is  nearly  identical  to  the  hazard  when  uncertain  parameters 
are  fixed  to  their  mean  values,  then  the  seismic  safety  of  individual 
facilities  can  be  calculated  through  conventional,  one-curve  analysis 
algorithms,  ignoring  parameter  uncertainty. 

One  may  conclude  that  both  types  of  hazard  analysis  (with  and  with¬ 
out  mixing  frequentist  and  non-frequentist  uncertainty)  are  useful, 
but  in  different  circumstances.  When  faced  with  a  specific  problem,  one 
should  first  ask  whether  quantified  non-frequentist  uncertainty  about 
frequentist  hazard  would  be  of  use.  If  the  answer  is  no,  as  it  seems  to 
be  the  case  for  the  safety  evaluation  of  given  facilities,  then  standard 
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procedures  with  best  estimates  of  the  parameters  are  usually  adequate. 

If,  on  the  other  hand,  decisions  depend  on  non-frequentist  uncertainty 
about  hazard,  as  it  is  the  case  in  research  resource  allocation,  then 
quantification  of  this  uncertainty  is  also  needed. 

There  are  cases  of  the  first  type  (only  the  average  hazard  is 
relevant)  that  cannot  be  solved  by  simply  performing  one  hazard  analysis 
with  parameters  fixed  to  mean  values,  because  the  results  would  be  biased 
or  because  it  is  difficult  to  define  or  work  with  mean  parameter  values. 

For  example,  it  may  be  impractical  to  define  a  mean  random  arrival  process, 
given  several  alternative  representations  of  earthquake  times,  and  it  may  be 
impossible  to  combine  the  opinions  of  different  experts  to  obtain  a  "mean 
opinion."  In  this  case,  expectation  of  hazard  should  be  taken  after  each 
model  has  been  analyzed  (possibly  in  combination  with  input  from  each  of 
several  experts).  Associated  with  each  model-expert  combination  there 
is  a  probability  and  a  conditional  hazard  curve  and  the  average  hazard  is 
obtained  by  weighing  the  conditional  hazard  functions  by  the  associated 
probabilities. 

Sensitivity  of  the  average  hazard  to  modeling  alternatives  depends  on 
the  site  under  consideration  and  on  the  expert  whose  opinion  is  being  con¬ 
sidered.  This  dependence  has  been  amply  documented  in  a  recent  study  by 
TERA  Corp.  (1980)  of  seismic  hazard  at  sites  of  the  central  and  eastern 
United  States.  The  same  study  shows  that  sensitivity  to  reasonable  param¬ 
eter  variations  depends  also  on  the  level  of  hazard  (the  return  period)  be¬ 
ing  considered  and  on  the  measure  of  site  intensity;  e.g.,  it  is  different 
for  peak  ground  acceleration  (PGA),  peak  ground  velocity  (PGV),  and  for  peak 
spectral  response  depending  on  the  natural  frequency.  It  is  therefore 
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clear  that  only  general  qualitative  conclusions  can  be  drawn  about  the 
influence  of  parameter  variations  and  errors  on  seismic  hazard.  When 
dealing  with  specific  situations,  these  conclusions  should  be  taken  as 
guidelines,  not  as  absolute  facts. 

Some  results  from  the  TERA  Corp.  (1980)  study  of  seismic  hazard 
in  the  central  and  eastern  United  States  will  be  discussed.  Nine  nuclear 
central  power  plant  sites  were  analyzed  using  input  information  from  ten 
different  experts.  Each  expert  responded  to  a  comprehensive  questionnaire 
on  seismic  source  zonation,  seismicity  of  each  identified  source,  and 
spatial  attenuation  of  motion  intensity,  provdinq  not  only  best  estimates 
of  the  parameters  but  also  quantified  uncertainty  on  i anges  of  values.  In 
the  analysis,  zonation  and  seismicity  parameters  (mean  rate  of  occurrence, 
decay  of  the  magni tude-recurrence  law,  and  upper  bound  magnitude)  were  allowed 
to  vary  from  expert  to  expert.  Three  different  attenuation  models  were 
used  and  kept  the  same  for  all  experts.  Two  of  them  have  the  form 

en(GM)  =  C1  +  C2I0  +  C3r  +  r  (4.2) 

in  which  GM  is  ground  motion  intensity  at  the  site  (PGA,  PGV,  or  spectral 
response),  the  coefficients  have  the  values  in  Table  4.1  for  the  at¬ 
tenuation  of  peak  acceleration  and  velocity,  and  r  is  epicentral  distance 
in  kilometers.  The  two  sets  of  coefficients  in  the  table  correspond  to 
using  different  MM  intensity  attenuation  relationships:  one, 

I$  =  IQ  +  3.2  -  O.OOllr  -  2.7  log1Qr  ,  I$  >  IV 

r  >  20  km 


(4.3) 
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Intensity 

Attenuation 

GM 

Units 

C1 

C2 

C3 

C4 

PGA 

cm/s^ 

1.18 

0.56 

-0.007 

-0.189 

Ossippee 

PGV 

cm/s 

-1.84 

0.64 

-0.006 

-0.167 

PGA 

cm/s2 

2.70 

0.56 

-0.0007 

-0.733 

Gupta-Nuttli 

PGV 

cm/s 

-0-50 

0.64 

-0.0006 

-0.645 

TABLE  4.1  Constants  in  the  attenuation  function 
of  Eq.  4.2  for  GM  =  PGA,  PGV  (after 
TERA,  1980). 
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was  adapted  from  Gupta  and  Nuttli  (1976);  the  other, 

Is  +  IQ  +  0.774  -  0.00117r  -  0.302  £n  r  (4.4) 

was  obtained  from  analysis  of  the  Ossippee  1940  earthquake  data.  A 
comparison  of  median  attenuation  models  is  shown  in  Figs.  4.1  and  4.2. 

The  figures  include  three  other  attenuation  relationships,  one  of  which 
(denoted  August  1979)  was  used  by  TERA  in  a  previous  evaluation  of  seismic 
hazard.  The  comparison  is  in  terms  of  sustained  acceleration  (=  0.7  PGA) 
for  magnitudes  5.5  and  6.5  and  uses  the  relationship  IQ  =  2mb  -  3.5, 
where  mb  =  body-wave  magnitude.  The  two  more  recent  models  (denoted  by 
"JAN  80"  in  the  figures)  intersect  around  20  and  200  km;  the  Ossippee 
attenuation  is  more  conservative  in  this  intermediate  range  of  epicentral 
distances  and  less  conservative  for  r  >  200  km. 

For  each  (median)  attenuation  function,  attenuation  uncertainty  was 
expressed  by  modeling  tn(GM)  as  a  truncated  normal  random  variable  with 
standard  deviation  a.  The  value  of  o  was  considered  to  be  either  0.7  or 
0.9  (values  around  0.6  result  from  direct  regression  of  acceleration  data 
in  the  west,  whereas  values  about  0.9  are  often  obtained  from  the  mathe¬ 
matical  combination  of  various  sources  of  uncertainty).  Use  of  the 
untruncated  normal  distribution  for  £n(GM)  is  not  satisfactory,  especially 
for  large  median  values  of  GM  and  large  o;  for  example,  in  the  case  when 
o  =  0.9,  a  3-o  departure  from  the  mean  of  £n(GM)  corresponds  to  an  un¬ 
realistically  high  factor  of  about  15  on  GM.  Truncation  of  the  distribu¬ 
tion  of  £n(GM)  was  therefore  introduced,  at  2  or  3  standard  deviations  from 
the  mean.  A  model  with  variable  a  was  also  considered,  with  a  value  of  0.6 
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FIGURE  4.1  Comparison  of  attenuation  models  for 
magnitude  5.5  (TERA,  1980) 
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for  r  <40  km  and  0.9  for  larger  distances.  Unfortunately,  there  is  little 
physical  or  statistical  basis  for  selecting  the  truncation  point,  which  is 
often  a  critical  parameter. 

Sensitivity  analyses  were  conducted  with  respect  to  possible  inter¬ 
pretations  of  each  expert  input  and,  in  all  cases,  with  respect  to  the 
median  attenuation  model  and  the  attenuation  error  distribution  mentioned 
above.  The  main  results  are  as  follows. 

In  nearly  all  cases,  sensitivity  to  upper  bound  magnitude  and  zonation 
alternatives  as  specified  by  each  expert  was  found  to  be  small.  Sensitivity 
to  the  upper  bound  magnitude  is  a  function  of  the  slope  of  the  magnitude- 
frequency  relationship  and  of  the  return  period.  If  the  return  period  is 
short  (e.y.,  1000  years  or  less)  and  the  slope  is  steep,  then  is 
usually  unimportant.  For  all  sites  located  in  the  Central  Stable  Region, 
sensitivity  to  zonation  alternatives  was  found  to  be  negligible  because 
the  main  contributing  sources  are  the  host  region  with  rather  well-defined 
boundaries,  and  the  New  Madrid  region,  which  is  too  distant  to  reflect 
minor  variations  in  its  boundaries.  Higher  sensitivity  was  noticed  at 
northeastern  sites  where  zonation  is  more  complex  and  zonation  uncertainty 
is  larger.  However,  also  in  this  case,  the  variation  of  hazard  when  using 
zonations  from  different  experts  exceeds  the  variation  associated  with 
zonation  uncertainty  of  individual  experts. 

The  median  attenuation  model  and  the  distribution  of  attenuation  error 
were  found  to  be  the  most  critical  parameters,  especially  for  PGV  and  long- 
period  spectral  ordinates.  Sensitivity  of  PGA  and  PGV  to  the  median  attenua¬ 
tion  at  four  sites  in  the  central  United  States  and  five  sites  in  the  eastern 
United  States  are  summarized  in  Tables  4.2  and  4.3,  respectively. 
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Sensitivity  varies  from  site  to  site  depending  on  its  relative  location 
with  respect  to  the  main  sources;  it  is  highest  for  P6V  and  for  sites 
where  hazard  is  largely  influenced  by  distance  sources. 

At  sites  in  the  east,  sensitivity  of  PGA  to  the  median  attenuation  is 
smaller  than  in  the  central  United  States  because  of  absence  of  important 
distant  sources  (the  St.  Lawrence  source  is  too  far  to  give  appreciable 
contributions).  However,  sensitivity  of  PGV  remains  high. 

The  distribution  of  the  attenuation  error  is  also  critical,  especially 
for  high  intensities  (long  return  periods).  For  a  return  period  of  1000 
years,  changing  the  truncation  point  from  two  to  three  standard  deviations 
increases  the  PGA  from  20  to  35  percent  if  o  =  0.9  and  from  15  to  25 
percent  if  o  =  0.7.  For  PGV,  these  quantities  are  respectively  8  to  15  per¬ 
cent  and  6  to  12  percent.  Higher  sensitivity,  up  to  45  percent,  is  found  for 
a  return  period  of  4000  years.  Reducing  the  attenuation  error  variance 
at  short  distances  (using  o  =  0.6  for  r  <  40  km  and  o  =  0.9  for  r  >  40  km, 
as  opposed  to  a  =  0.9  for  all  r)  produces  very  minor  changes  in  the  hazard, 
except  for  the  few  sites  where  hazard  depends  on  near  sources. 

As  one  would  expect,  the  importance  of  far-away  sources  increases 
with  increasing  period  (from  PGA  to  PGV)  and  the  importance  of  sources  with 
higher  upper  bound  magnitude  increases  with  decreasing  hazard.  Also,  the 
relative  importance  of  different  sources  depends,  in  some  cases  rather 
importantly,  on  the  form  of  the  attenuation  function. 

Average  hazard  and  sensitivity  results  for  firm-soil  sites  in  Boston 
have  been  obtained  by  Cornell  and  Merz  (1975).  The  zonation  for  the  base 
case  is  shown  in  Fig.  4.3  and  best  estimates  of  the  seismicity  parameters 
for  each  source  are  given  in  Table  4.4.  The  attenuation  law,  in  terms  of 
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SOURCE 

No.  of  events 
with 

MM I  >  4.25/year/10,000  km? 

Slope 

Upper  bound  intensity 

1 

3.16 

0.48 

7.7 

2 

0.42 

0.48 

8.7 

3 

0.180 

0.48 

6.7 

4 

0.150 

7.3 

5 

0.010 

0.48 

6.3 

6 

0.215 

0.48 

6.7 

7 

0.200 

0.48 

6.7 

8 

0.380 

0.48 

7.3 

TABLE  4.4  Best  estimates  of  source  parameters  for  hazard 
analysis  of  Boston  (Cornell  and  Merz,  1975). 


Modified  Mercalli  intensity,  was  estimated  from  the  isoseismals  of  large 
historical  earthquakes  in  the  region;  it  has  the  form 


(In  ,  r  <  10  miles 

°  (4.5) 

2.6  +  Ig  -  1.3£nr,  r  >  10  miles 

Attenuation  intensity  is  expressed  by  an  additive  normal  random  variable 
with  mean  0  and  standard  deviation  0.2.  The  main  sensitivity  cases 
considered  by  Cornell  and  Merz  correspond  to  the  following  variants  of 
the  base  case: 

1.  Is  increased  or  decreased  by  1/2  unit; 

2.  Source  8  replaced  with  a  more  conservative  or  a  less 
conservative  source. 

3.  The  upper  bound  epi central  intensity  of  all  sources  increased 
or  decreased  by  1/2  unit. 

4.  The  upper  bound  epicentral  intensity  for  sources  1  and  2  (off 
Cape  Ann)  increased  to  XII. 

Hazard  curves  corresponding  to  these  and  other  parameter  perturbations 
are  shown  in  Fig.  3.1.  Decreasing  site  intensity  by  1/2  unit  in  the 
attenuation  equation  (curve  AL)  produces  appreciable  changes  only  for 
short  return  periods.  For  long  return  periods  the  hazard  remains  the 
same  because  the  contributing  events  are  primarily  those  very  close  to 
the  site,  within  a  10-mile  radius.  Conservative  modification  of  Source  8 
(curve  GP)  produces  similar  qualitative  results,  with  noticeable  effects 
only  for  short  return  periods.  More  important  sensitivity  is  found  with 
respect  to  the  upper  bound  intensities.  If  these  bounds  are  increased 
by  1/2  unit  (curve  UBP),  then  the  hazard  curve  moves  to  the  right  by  an 
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amount  of  up  to  0.4  for  long  return  periods.  Larger  increases  of  intensity 
over  long  return  periods  are  observed  if  the  Cape  Ann  sources  are  assigned 
upper  bound  intensity  XII  {curve  UB12). 

For  calculation  of  the  composite  (average)  hazard  curve,  probabilities 
were  subjectively  assigned  to  discrete  alternatives  about  the  parameters  of 
the  attenuation  function,  the  geometry  of  Source  8,  and  the  upper  bound 
intensities.  These  probabilities  are  given  in  Table  4.5.  They  correspond 
to  36  different  cases.  The  resulting  composite  hazard  curve  is  compared 
in  Fig.  4.4  with  the  curve  for  parameters  fixed  to  most  likely  values. 

The  very  significant  divergence  for  return  periods  larger  than  10,000 
years  is  due  to  the  possibility  that  the  maximum  intensity  is  XII  for 
the  sources  near  Cape  Ann. 

Additional  sensitivity  analyses  for  Boston  have  been  conducted  by 
TERA  (1980)  with  respect  to  expert  input,  median  attenuation  model  (Eq. 

4.4  vs.  Eq.  4.5),  and  dispersion  of  the  attenuation  error  (o  =  1.0  in 
place  of  o  =  0.2).  The  results  are  summarized  in  Table  4.6  for  return 
periods  of  1,000  and  10,000  years.  As  shown  in  Fig.  4.5,  the  Gupta- 
Nuttli  attenuation  model  produces  site- intensity  values  about  1/2  unit 
higher  than  the  Cornell-Merz  attenuation.  The  effect  on  hazard  is  a 
variation  of  about  0.3  units  of  intensity.  More  sensitive  parameters 
are  the  expert  opinion  and  the  dispersion  of  the  attenuation  error. 

Although  the  expert  in  the  TERA  analysis  basically  agreed  with  the 
zonation  by  Cornell  and  Merz,  he  assigned  different  values  to  upper 
bound  intensities  and  occurrence  rates,  with  a  net  effect  of  decreasing 
the  hazard  by  approximately  0.9  units  of  intensity.  An  even  more 
critical  choice  is  that  of  o  which,  when  varied  from  0.2  to  1.0, 
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LOW 

BEST 

ESTIMATE 

HIGH 

MM I  XII 
for  Cape  Ann 

+ 

ATTENUATION  LAW 

0.25 

0.50 

0.25 

ZONATION 

0.25 

0.50 

0.25 

MM I  UPPER  BOUNDS 

0.20 

0.30 

TABLE  4.5  Subjective  probabilities  for  the  calculation  of 

the  composite  hazard  curve  (Cornell  and  Merz,  1975). 
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FIGURE  4.4  Most  likely  and  composite  hazard  curves  for 
Boston  (Cornell  and  Merz,  1975). 
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produces  an  increase  of  site  intensity  of  1.1  and  1.45  units,  for  return 
periods  of  1,000  and  10,000  years,  respectively. 

Sensitivity  results  for  the  La  Crosse  site  (Fig.  4.6)  from  Dames  and 
Moore  (1980)  and  TERA  (1980)  are  given  next.  Sensitivity  involves 
zonation,  recurrence  parameters,  and  median  attenuation.  In  all  cases, 
the  attenuation  error  on  £n(PGA)  is  modeled  as  an  untruncated  normal 
variable  with  standard  deviation  a  =  0.6.  Two  alternative  zonations 
are  considered:  the  Central  Stable  Region  only,  or  the  Central  Stable 
Region  and  the  Wisconsin  Arch  source  shown  in  Fig.  4.6.  In  each  case, 

sensitivity  is  evaluated  also  with  respect  to  inclusion  or  exclusion  of 

the  New  Madrid  source.  The  median  attenuation  functions,  based  on  the 
Ossippee  and  Gupta-Nuttli  models,  are  the  same  as  described  earlier  in 
this  section.  Table  4.7  gives  the  seismicity  parameters,  as  assigned  by 
Dames  &  Moore  (D&M)  and  by  TERA  (Eg).  Especially  noticeable  in  the  table 
is  the  difference  in  mean  activity  rate  (which  is  much  higher  in  the  TERA 
report)  and  in  the  upper  bound  magnitude  (the  three  values  in  the  last 

column  correspond  to  the  best  estimate  and  to  lower  and  upper  bounds). 

The  difference  in  the  best  estimate  is  about  1/2  magnitude  unit  for  all 
sources. 

Hazard  results  in  Table  4.8  show  that  the  effects  on  PGA  of  changing 
the  mean  rate  a  and  the  upper  bound  magnitude  M^  are  about  15  and  10 
percent,  respectively.  Sensitivity  to  the  inclusion  or  exclusion  of  the 
New  Madrid  source  is  entirely  dependent  on  the  attenuation  model:  it  is 
negligible  for  Ossippee  and  about  25  percent  for  Gupta-Nuttli.  Finally, 
the  effect  of  including  or  excluding  the  Wisconsin  Arch  source  is  not 
large  but  about  the  same  for  both  attenuation  functions. 
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FIGURE  4.6  La  Crosse  site  and  the  Wisconsin  Arch 
source  (Dames  and  Moore,  1980). 


Dames  and  Moore  (1980)  and  TERA  (1980). 
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PGA  (g)  AT  LA  CROSSE  FOR 
1000  AND  10,000  YEAR  RETURN  PERIOD 


Ossippee 

Gupfo-Nuttli 


Central  Stable  Region 
without  New  Modrid 

Central  Stoble  Region 
and  New  Modrid 

DAM  a  value 

DAM  M 

D&M  a  value 
D&M  My 

E.  a  value 
DAM  My 

Eq  a  value 
E>u 

• 

E_  a  value 

E9Mu 

RP  1 

RP  2 

RP  1  RP  2 

RP  1  RP  2 

RP  1  RP  2 

RP  1  RP  2 

0.061 

0.15 

0.061  0.15 

0.081  0.18 

0.085  0.20 

0.085  0.20 

0.052  0.12 

0.065  0.14 

0.075  0.15 

0.076  0.17 

0.079  0.17 

New  Madrid 
Excluded 

Centrol  Stable  Region  and  Northern  Illinois  (Eq) 
and  New  Madrid 

D&M  o  value 
D&M  M 

u 

D&M  a  value 
D&M  M 

u 

E„  o  volue 
Duivi  M^ 

Eq  a  value 

Vu 

Eq  o  value 

RP  1 

RP  2 

RP  1  RP  2 

RP  1  RP  2 

RP  1  RP  2 

RP  1  RP  2 

0.058  0.12 

0.058  0.12 

0.068  0.14 

0.084  0.17 

0.085  C.I7 

0.049  0.096 

0.062  0.1? 

0.066  0.13 

0.074  0.15 

0.075  0.15 

Ossippee 

Gupfo-Nuttli 


Considering  mb=4.0  and  greater 

RP  1  =  1,000  year  return  period 
RP  2  =  10,000  year  return  period 


TABLE  4.8  Sensitivity  results  for 

La  Crosse  site  (TERA,  1980) 


77 


All  parameter  sensitivity  results  up  to  here  refer  to  essentially 
identical  seismic  hazard  models.  If  the  model  is  changed,  then  not  only 
the  numerical  value  of  estimated  hazard  will  change  (see  Sec.  4.3)  but 
also  sensitivity  to  parameter  variations  may  be  different.  One  example 
is  the  line  source  model  of  Der  Kiureghian  and  Ang  (1975).  The  main 
departure  from  the  standard  model  is  that  earthquakes  are  now  associated 
with  fault  ruptures  of  finite  dimension  and  site  effects  are  related  to 
the  (minimum)  distance  to  the  rupture  zone  as  opposed  to  the  epicentral 
or  focal  distance.  As  a  consequence  of  this  modification,  calculated 
hazard  increases  and  parameter  sensitivity  changes.  With  regard  to  the 
latter  effect,  Der  Kiureghian  and  Ang  challenge  the  conclusions  by 
Cornell  and  Vanmarcke  (1969)  based  on  point-source  models,  that 

1.  The  major  contribution  to  hazard  comes  from  frequent,  small -size 
earthquakes  from  sources  close  to  the  site. 

2.  Hazard  is  seldom  sensitive  to  the  upper  tail  of  the  magnitude 
distribution  or  to  the  upper  bound  magnitude. 

From  the  line-source  model,  Der  Kiureghian  and  Ang  find  that 

1.  Distant  sources  may  have  significant  effects  on  seismic  hazard 
when  large-magnitude  earthquakes  are  possible. 

2.  Seismic  hazard  is  sensitive  to  the  upper  bound  magnitude  My, 
particularly  for  My  <  8.5. 

These  conclusions  are  especially  significant  for  long  return  periods. 

The  main  reason  why  point  and  line  source  models  have  different  sensitivity 
to  high  magnitude  earthquakes  is  that  these  events  are  usually  associated 
with  long  rupture  zones,  which  may  propagate  close  to  the  site  and  produce 
effects  much  larger  than  expected  from  epicentral  distance. 
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For  practical  applications,  Der  Ki'ureghian  and  Ang  conclude  that 
the  choice  of  My  is  inmaterial  at  low  intensity  levels,  that  in  all  cases 
events  of  extremely  large  size  (magnitude  larger  than  9)  give  negligible 
contribution  to  hazard,  and  that  the  choice  of  My  is  critical  only  for 
long  return  periods  and  when  contributing  sources  exist  far  from  the 
site  with  My  below  8.5.  They  also  conclude,  somewhat  contrary  to  findings 
in  the  TERA  study,  that  sensitivity  to  My  is  not  much  affected  by  the  value 
of  the  slope  parameter  b. 

Fig.  4.7  shows  a  comparison  of  the  point  and  line  source  models  and 
some  sensitivity  results  for  the  latter  model.  All  curves  in  the  figure 
estimate  hazard  in  San  Francisco  and  neglect  attenuation  uncertainty. 

With  reference  to  the  line-source  model,  curve  1  represents  the  base  case 
evaluation;  curve  2  results  from  neglecting  events  with  magnitude  below 
4.5;  and  curve  3  corresponds  to  increasing  the  slope  of  the  magnitude 
recurrence  law  from  1.22  to  1.30.  The  remaining  two  curves  account  for 
incompleteness  of  the  data  when  fitting  the  magnitude  recurrence  relation¬ 
ship  (curve  4)  and  reinterpret  the  rate  of  activity  along  the  San  Andreas 
fault  (curve  5).  The  last  curve  may  be  considered  to  quantify  sensitivity 
to  zonation.  In  all  cases,  the  upper  bound  magnitude  was  fixed  to  8.5. 
Sensitivity  is  clearly  small,  especially  when  compared  with  the  effect 
of  including  or  excluding  attenuation  uncertainty.  The  latter  effect  and 
the  importance  of  the  choice  of  the  median  attenuation  curve  is  demonstrated 
in  Fig.  4.8  for  which  the  three  attenuation  models  of  Table  4.9  have  been 
used,  with  and  without  uncertainty.  It  is  interesting  to  note  that 
models  2  and  3,  which  are  substantially  different  ir  their  median  pre¬ 
dictions,  produce  nearly  identical  results  when  attenuation  uncertainty 
is  included,  whereas  the  opposite  is  true  for  models  1  and  2. 
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FIGURE  4.7  Comparison  of  point  and  line  source  models 
and  sensitivity  analysis  for  San  Francisco 
(Der  Kiureghian  and  Ang,  1975). 
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FIGURE  4.8  Sensitivity  of  the  line-source  model  to  median  attenuation 
and  to  attenuation  uncertainty.  Hazard  estimates  are  for 
San  Francisco  (Der  Kiureghian  and  Ang,  1975). 
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Eq. 

MEDIAN  ATTENUATION 

a 

1 

a  =  1320e0,58m  (r+25)'1,52 

0.894 

2 

a  =  I080e°-50m  (r+25)'1,32 

0.755 

3 

a  =  472.3e0,64,n  (r+25) _1  * 301 

0.573 

TABLE  4.9  Attenuation  models  used  in  the 

sensitivity  analysis  of  Fig.  4.8. 
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The  exact  determination  of  source  boundaries  is  usually  not  of  critical 
importance  for  "area"  sources,  but  there  are  notable  exceptions.  One  ex¬ 
ception  is  when  different  interpretations  are  possible  for  the  mechanism  of 
earthquake  generation.  For  example,  it  has  been  suggested  by  Long  (1976) 
and  McKeown  (1978)  that  high  intensity  seismicity  in  some  regions  of  the  east 
ern  United  States  might  be  due  to  (or  might  be  statistically  associated  with) 
the  presence  of  plutons.  More  frequent  situations  in  which  source  geometry 
plays  an  important  role  are  when  it  is  uncertain  whether  a  major  fault  is 
active  or  not.  A  study  summarized  in  Donovan  and  Bornstein  (1978)  shows 
that  the  change  of  hypothesis  about  activity  or  inactivity  of  a  fault  in  the 
Los  Angeles  areas  leads  to  error  factors  of  up  to  50  percent  on  site  accelera 
tion  for  a  given  exceedance  probability.  This  should  be  regarded  as  a 
rather  extreme  case.  Sensitivity  of  hazard  to  upper  bound  magnitude  depends 
on  the  value  of  the  maximum  magnitude  and  on  the  location  of  the  site 
relative  to  the  source,  but  it  is  now  recognized  that  this  parameter  is  not 
as  critical  as  previously  believed  (see,  for  example,  Cornell,  1980;  McGuire 
and  Shedlock  ,  1980;  Donovan  and  Bornstein,  1978). 

The  diversity  of  conclusions  from  sensitivity  studies  in  this  section 
demonstrates  the  impossibility  of  giving  absolute  rules  about  the  influence 
of  parameter  uncertainty  on  hazard.  At  most,  one  can  identify  sources  of 
uncertainty  that  tend  to  dominate  the  results;  among  these  sources  are  the 
form  of  the  attenuation  function  and  the  distribution  of  the  attenuation 
error.  Expert  opinion  should  also  be  included  as  a  major  source  of  hazard 
variation,  as  will  be  said  in  Sec.  6. 
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4.3  Alternative  Models 

Alongside  the  development  of  what  we  now  regard  as  standard  seismic 
hazard  procedures,  many  alternative  modeling  hypotheses  have  been  explored. 

Much  work  has  been  done  on  the  representation  of  earthquake  times 
(sometimes  also  locations)  through  random  point  processes  other  than  Poisson. 
The  motivating  idea  is  that  earthquakes  produce  changes  in  the  state  of 
the  physical  system  and  must  therefore  affect  the  occurrence  of  other 
seismic  events  at  close  times  and  locations.  Perhaps  the  simplest  way  to 
introduce  dependence  of  earthquakes  in  time  is  through  so-called  renewal 
and  clustering  models. 

In  a  renewal  point  process,  interarrival  times  T-| ,  T2,  ...  are  indepen¬ 
dent  and  identical  variables  as  in  the  Poisson  model,  but  not  necessarily 
with  exponential  distribution.  Therefore,  memory  extends  only  between 
consecutive  events.  Kameda  and  Ozaki  (1979)  fitted  a  simple  renewal  model 
to  historical  data  from  Kyoto.  They  found  that  these  data  are  compatible 
with  a  piecewise  constant  hazard  function  v(t)  =  fj(t)/[l-Fj(t)],  with 
only  two  values:  Vj  for  small  t  and  for  large  t, and  compared 

hazard  predictions  from  this  model  with  those  when  v(t)  =v (Poisson  events). 
A^ter  observing  that  in  the  Kyoto  catalogue,  longer  interarrival  times 
usually  terminate  in  larger  earthquakes,  Kameda  and  Ozaki  generalized  the 
model  to  account,  for  this  dependence  and  again  compared  results  with 
the  standard  assumptions  of  Poisson  arrivals  and  independent  sizes.  The 
difference  in  hazard  is  sometimes  large,  depending  on  the  interval  of 
prediction  and  on  the  time  of  occurrence  of  the  last  earthquake.  One  set 
of  results  is  shown  in  Fig.  4.9.  Each  curve  in  the  finure  gives  the 
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FIGURE  4.9  Probability  of  future  earthquakes  of 
JMA  intensity  VI  or  greater  in  Kyoto, 
starting  from  year  1978.  Poisson 
and  renewal  models.  (Kameda  and 
Ozaki,  1979). 
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probability  of  at  least  one  earthquake  in  Kyoto  with  JMA  (Japan  Meterological 
Agency)  intensity  VI  or  greater  (approximately  MM  intensity  IX  or  greater) 
as  a  function  of  time,  starting  from  the  year  1978.  The  dashed-dotted 
line  corresponds  to  the  hypothesis  of  Poisson  events  with  independent  sizes. 
The  other  two  lines  are  for  arrivals  according  to  the  renewal  process:  the 
dashed  line  continues  to  assume  independent  sizes,  whereas  the  solid  line 
incorporates  dependence  of  size  on  interarrival  time.  In  the  latter  case, 
the  distribution  of  intensity  has  either  of  two  forms,  depending  on  whether 
the  earthquake  occurs  when  v  =  v-|  or  v  =  xu.  It  is  interesting  to  observe 
that  in  the  present  case  and  for  short  to  moderate  time  intervals,  dependence 
of  size  on  interarrival  time  counteracts  the  effect  on  hazard  of  non-Poisson 
arrivals.  A  different  renewal  process  of  earthquake  times  was  proposed 
earlier  by  Tsteva  (1976). 

Clustering  models  have  been  developed  more  extensively  (see  for  example, 
Vere-Jones  and  Davies,  1966;  Walley,  1976;  Schlien  and  Toksoz,  1970;  Knopoff 
and  Kagan,  1977).  Each  cluster  comprises  one  "main  shock"  and  the  associ¬ 
ated  family  of  "aftershocks."  Clusters  are  typically  assumed  to  occur  at 
Poisson  times  (and  Poisson  locations  in  space),  in  spite  of  some  evidence 
of  secondary  "clustering  of  clusters"  over  long  times  and  distances  (Walley, 
1976).  Clustering  has  of  course  the  effect  of  increasing  hazard  at  sites 
and  times  close  to  sites  and  times  of  past  earthquake  occurrences.  However, 
it  is  found  that  if  the  main  event  is  also  the  most  intense  event  of  the 
cluster  and  aftershocks  are  confined  to  a  small  region  near  the  epicenter 
of  the  parent  earthquake,  then  the  hazard  consequence  of  excluding  after¬ 
shocks  is  small  (Merz  and  Cornell,  1973). 

An  interesting  extension  of  the  concept  of  uncertain  seismicity  rate 
A  has  bean  explored  by  Walley  (1976),  who  considers  as  random  the  intensity 
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function  A(X-|,  X ^ .  t)  of  the  Poisson  process  of  epicentral  coordinates 
X-j  and  Xg  and  times  t.  This  assumption  transforms  the  Poisson  process 
into  a  so-called  doubly-stochastic  Poisson  process.  Realizations  of  doubly- 
stochastic  Poisson  processes  are  very  similar  to  realizations  of  clustering 
models.  The  effect  of  seismic  hazard  is  also  similar.  In  an  application 
to  New  Zealand,  the  same  author  finds  that  hazard  is  sensitive  to  the  depth 
of  energy  release  (a  similar  conclusion  for  sites  near  major  faults  is 
reach**  by  McGuire  and  Shedlock  ,  1980)  and  cautions  that  significant  errors 
may  come  from  the  assumption  of  stationary  seismicity. 

On  the  basis  that  major  earthquakes  release  large  amounts  of  strain 
energy  and  that  such  energy  gradually  accumulates  in  time,  one  would  expect 
that  the  interarrival  time  be  related  to  the  magnitude  of  the  previous  and 
the  next  earthquake.  Convenient  mathematical  representations  of  this  depen¬ 
dence  have  been  given  in  terms  of  Markov  and  semi -Markov  processes  (Patwardhan 
et  al_. ,  1978,  Kiremidjiar.  et  a/L ,  1980).  Hazard  predictions  are  sometimes 
very  different  from  those  using  the  Poisson  model,  depending  on  the  time  and 
magnitude  of  the  last  major  earthquake. 

As  remarked  previously  in  Sec.  4.2,  hazard  is  usually  sensitive  to  the 
attenuation  function,  and  even  more  to  the  variance  and  upper  tail  distri¬ 
bution  of  the  attenuation  error.  The  variance  depends  not  only  on  the  data 
set  and  on  which  quantities  are  used  to  characterize  source  and  site  intensity, 
but  also  on  the  method  by  which  the  data  are  analyzed.  Multistep  procedures 
which  utilize  prediction  of  intermediate  parameters  lead  to  particularly 
large  estimates  of  dispersion  (Cornell  et  ah ,  1979). 

Errors  on  hazard  at  New  Zealand  sites  from  neglecting  azimuth-dependence 
of  the  attenuation  are  reported  in  Walley  (1976).  He  finds  that  when  an 
isotropic  model  is  replaced  with  a  set  of  highly  eccentric  elliptical  iso- 
seismals,  hazard  at  some  sites  varies  by  0.6  to  0.8  units  of  MM  intensity. 
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Of  course,  the  effect  of  anisotropy  is  less  pronounced  if  seismic  activity 
is  uniform  over  large  regions. 

Sensitivity  of  hazard  to  including  or  neglecting  the  spatial  dimension 
of  fault  rupture  (point  versus  line  source  model)  has  already  been  reported 
in  Sec.  4.2;  see  for  example  Fig.  4.7. 
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5.  ERRORS  FROM  MODELING  ASSUMPTIONS 

As  said  earlier  in  Sec.  1.2,  the  correct  but  unknown  model  MT  does 
not  necessarily  belong  to  the  family  of  models  ®.  The  latter  comprises  only 
simple  models  to  which,  for  convenience  of  analysis,  My  is  assumed  to  belong. 

For  example,  all  models  in®  might  assume  isotropic  attenuation  on  the  geo¬ 
graphical  plane,  whereas  actual  attenuation  may  be  azimuth-dependent  due  to 
oriented  geological  structures  or  to  the  directional  propagation  of  fracture. 

One  should  expect  that  if  My  bias  will  result  from  assuming  the  contrary. 

If  bias  is  small,  then  the  seismic  hazard  procedure  is  said  to  be  "robust." 

Robustness  is  a  desirable  property  because  My  is  unknown  and  certainly  not 
of  the  elementary  type  hypothesized  for  analysis.  In  what  follows,  a  true 
model  My  is  assumed  and  a  few  common  estimators  of  seismic  hazard,  H,  are 
evaluated  against  the  associated  true  hazard,  Hy.  Special  attention  is 

given  to  the  conditional  bias  (HT  -  E  [H])  in  which  Z  is  statistical 

T  ZlMj 

sampling  information.  In  reality,  the  analyst  is  ignorant  of  My  and  Hy, 
his  only  statistical  information  being  a  historical  sample  Z,  e.g.,  of  ground 
accelerations  at  the  site.  The  advantage  of  knowing  Hy  is  that  it  will  be 
possible  here  to  evaluate  bias  and  error  variance  of  various  seismic  hazard 
procedures  which  the  analyst  might  choose  to  use.  Contrary  to  findings  by 
McGuire  and  Shedlock  (1980)  for  the  case  when  My  belongs  to®  (see  Sec.  4.1), 
error  bias, not  variance, is  found  to  be  the  primary  source  of  mean  square 
error  of  the  hazard  estimators,  especially  at  high  levels  of  intensity. 

Material  in  this  section  is  for  the  most  part  from  Schumacher  (1977)  and 
from  unpublished  work  by  Veneziano  et  al .  (1977). 

5.1  True  Model  and  Hazard  Estimators 

Suppose  that  the  true  hazard  at  the  site,  Hy,  corresponds  to  the  fol¬ 
lowing  probabilistic  representation  of  regional  seismicity,  My: 

I 
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a.  There  is  only  one  source,  in  the  form  of  a  disc  with  the  center  at  the 

site  and  radius  r.  Earthquakes  occur  as  Poisson  points  in  time,  with  a 

2 

mean  annual  rate  of  vt  =  virr  while  in  space;  epicenters  are  uniformly 
and  independently  distributed  inside  the  source. 

b.  Earthquake  magnitudes  Mp  M2,  ...  are  independent  random  variables 
with  identical,  doubly-truncated  exponential  distribution 

-b(m-mQ) 

FM(m)  =  l<i  [1  -  e  ]  .  i"q  <  m  $  m-j  (5.1) 

in  which  mQ  and  m-j  are  lower  and  upper  bounds  to  M,  b  is  a  positive 
constant,  and  K-j  =  {1  -  exp[-b(m.|-mQ]}"^  is  a  normalization  constant. 

c.  Peak  acceleration  at  the  site,  a,  is  the  only  intensity  parameter  of 
interest.  It  is  related  to  magnitude  M  and  epicentral  distance  R  as 
(Esteva,  1970;  Cornell,  1971;  Der  Kiureghian  and  Ang,  1977) 

fc-M  -b, 

a  =  b1  e  (R+b4)  *  c  (5.2) 

in  which  the  bi  are  constant  parameters  and  e  is  a  lognormal  random 
variable  with  independent  realizations  for  different  earthquakes. 

d.  Hazard  is  measured  in  terms  of  the  probability  G(a,  t)  that  ground 
acceleration  at  the  site  will  exceed  a, at  least  once  in  t  years. 

Under  these  assumptions,  the  correct  mathematical  representation  of  hazard, 
Hy  =  Gy(a,  t),  is  given  by 

-vftG  (a) 

HT  =  Gy(a,  t)  =  e  1  a  (5.3) 
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in  which  G„(a)  is  the  correct  complementary  CDF  of  peak  ground  acceleration 

a 

at  the  site  due  to  an  event  of  random  magnitude  and  epicentral  location. 

The  function  G_ (a )  has  not  simple  analytical  form,  but  it  can  be  calculated 

a 

numerically,  e.g.,  by  the  method  of  Cornell  and  Merz  (1975). 

Given  a  finite  record  of  peak  ground  accelerations  at  the  site,  the 
quantities  vt  and  Ga(a)  in  Eq.  5.3  cannot  be  calculated  exactly.  For  values 
of  a  of  practical  interest,  the  error  in  estimating  Gy(a,  t)  is  contributed 
primarily  by  errors  in  the  estimation  of  G_ (a ) .  Therefore,  it  will  be  assumed 

a 

here  that  v.  is  known  and  inference  will  be  limited  to  G,(a).  Moreover,  one 

U  fl 

may  limit  consideration  to  acceleration  values  above  a  threshold  a„.  of 

min 

engineering  interest  and  estimate  the  conditional  complementary  CDF, 

Ga  I  a  >  a  .  ^  =  Ga^  Ga^amin)  (5-4) 

'  min 

for  a  >  amin.  This  last  function  is  plotted  in  Fig.  5.1  (dashed  lines)  for 
the  following  parameter  values: 

'  amin  =  10  (cm/sec2) 
r  =  100  km 

mQ  =  4,  m.|  =  8,  b  =  2  (5.5) 

b1  =  463.2,  b2  =  0.64,  b3  =  1.301,  b4  =  25 
£n  e  ~N(0,  ole)  ,  o£nc  =  0,  0.6 

The  values  of  bp  b2,  b3,  and  are  from  Esteva  (1976).  It  is  further 
noticed  that  inference  of  the  complementary  CDF  of  a,  given  that  a  >  am-n  , 
is  equivalent  to  inference  of  the  complementary  CDF  of  any  of  the  following 
random  variables: 


log,n{Probability  of  exceedance)  given  an  event  with  a  a  10  cm/sec 
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FIGURE  5.1  True  hazard  curves  and  estimated  hazard  curves 
for  H  =  untruncated  exponential  family  and 
H  =  maximum-likelihood  estimator  of  Ga- 
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Y  =  An  a  (5.6) 

W  =  An  a  -  Jin  a  ..  ,  M  >  0 

It  is  convenient  here  to  characterize  estimators  of  hazard  as  inference 
rules  on  Gy(y)  or  on  G^(w). 

Estimators  of  Hy  =  G^(w)  may  differ  because  of  differences  in  the 
class  H  of  hazard  functions  to  which  Gj.,  is  assumed  to  belong  or  because  of 
different  inference  rules,  e.g,  maximum  likelinood,  least  squares,  Bayes. 
Both  aspects  have  important  effects  on  H  =  G^(w),  especially  for  large  w. 

In  selecting  a  specific  estimator  one  looks  for  simplicity,  richness  (H 
should  contain  a  variety  of  hazard  functions  so  that  at  least  one  is  close 
to  the  unknown  true  function  Gw),  and  accuracy.  The  difference  e(w)  = 

Gw(w)  -  Gw(w)  is  a  random  process  with  the  conditional  second  moment 
function 

Be(w)  =  E[e2(w)|Gw]  (5.7) 

This  function  should  be  as  small  as  possible,  at  least  for  a  certain  class 
of  true  complementary  CDF's,  G^. 

Richness  and  accuracy  are  related.  If  H  is  a  narrow  family,  then  Gy 
might  be  very  different  from  all  the  functions  inH.  As  a  consequence,  Gw 
is  likely  to  be  heavily  biased  and  Bg(w)  to  be  large.  At  the  other  extreme, 
if  H  is  very  rich,  then  bias  may  be  small  but  the  variance  of  e(w)  is  large 
and  Bg(w)  is  again  large.  A  good  choice  of  H  is  one  that  achieves  a  balance 
between  bias  and  error  variance  for  a  given  class  of  actual  hazard  functions 
and  one  that,  at  the  same  time,  is  easy  to  analyze.  Bias  and  error  variance 
have  been  calculated  for  the  following  estimation  rules  and  parametric 
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families  of  distributions: 

Maxi mum- likelihood  Estimation  wi th 
M  =  one-parameter  gamma  family. 

K  =  two-parameter  truncated-exponential  family. 

H  =  two-parameter  Wei  bull  family. 

K  =  two- parameter  lognormal  family. 

Bayesian  Estimation  with 

H  =  one-parameter  exponential  or  two-parameter  truncated-exponential 
family. 

R  =  two-parameter  Weibull  family. 

R  =  two-parameter  lognormal  family. 

In  some  other  cases,  hazard  at  the  site  was  estimated  indirectly  by  first 
fitting  a  regional  (as  opposed  to  site)  model  (source  configuration,  seis¬ 
micity  parameters,  attenuation  function)  to  statistical  data.  A  brief 
account  of  these  more  complicated  analyses  is  given  in  Sec.  5.4. 

For  all  site  models,  sample  information  was  assumed  in  the  form  of 
an  independent  sample  of  size  N,  2  =  (W-j ,  ...»  W^}  from  the  true  distribu¬ 
tion  of  W,  and  estimation  error  statistics  were  obtained  through  repeated 
Monte  Carlo  simulation  of  Z.  Specifically,  three  sets  of  independent  samples 
were  used:  Set  1  =  20  samples  each  of  length  N  =  20  from  the  distribution 
of  W  that  corresponds  to  o£ne  =  0;  Set  2  =  same  as  set  1  but  with  = 

0.6;  Set  3=17  samples  each  of  length  N  =  50  from  the  distribution  of  W 
that  corresponds  to  oAn  =  0.6.  Some  approximate  analytical  results  are 
also  given. 
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5.2  Performance  of  Simple  Maximum-Likelihood  Estimators 

Given  H  as  a  set  of  complementary  CDF's  and  given  a  sample  Z  from  Gy, 
the  maximum-likelihood  estimator  Gy  is  the  function  in  M  for  which  the 
probability  (density)  of  Z  is  maximum.  Detailed  analysis  follows  for  the 
parametric  families  H  given  previously  in  Sec.  5.1. 

a.  K  =  one-parameter  family  of  gamma  distributions 

Let  H  =  {GA(K,  X),  K  =  shape  parameter,  known}  be  a  parametric  family 
of  gainna  complementary  CDF's  with  X  as  the  only  parameter.  Consider  first 
the  case  when  K  =  1,  i.e.,  when  GA(K,  X)  is  the  exponential  function 

Gy(w)  =  e'Xw  ,  w  >  0  (5.8) 

Then,  as  a  function  of  Z  *  {W^,  ....  Wy),  the  maximun-likelihood  estimator 
of  X  is 


X  =  1/ff  ,  in  which  W  =  ^  £  W. 


(5.9) 


An  accurate  and  convenient  approximation  to  the  distribution  of  W  is  the 
gamma  distribution  with  (exact)  mean  m^  =  my  and  (exact)  variance  =  Oy/N. 
In  fact,  this  would  be  the  correct  distribution  of  if  W  itself  had  gainna 
distribution,  as  opposed  to  one  of  the  dashed  distributions  in  Fig.  5.2. 


Under  this  approximation,  -tnGy(w)  =  Xw  has  so-called  inverted  gamma-1 
distribution  (Raiffa  and  Schlalfer,  1961)  and  after  a  few  substitutions 
one  finds 


Peak  site  acceleration,  a (cm/sec  ) 


FIGURE  5.2  True  hazard  (dashed  line)  and  fractlles  of 
the  maximum-likelihood  estimator  (solid 
lines)  when  W  is  assumed  to  have  gamma 
distribution  with  shape  parameter  K  ■  3. 
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SC-AnGjw)]  = 


.  w  N 


"w  n-v; 


w 


.  N  >  V,' 


N2V2 


(5.10) 


Var[-“V“,]  =  %  • N  >  "5 


in  which  is  the  coefficient  of  variation  of  W.  The  mean  square  error  of 
the  maxi mum- likelihood  estimator  of  tnGy(w)  is 


Be(w)  =  bg(w)  +  Var[-S.n  Gw(w)] 


(5.11) 


N2V2 
N  VW 


"Yl  N-Vf 


m^  (N-V2)2(N-2V2) 


N  >  2V, 


W 


in  which  bg  denotes  bias.  Plots  of  the  expected  value  and  of  the  +1 
standard  deviation  interval  of  £n  Gw(w)  are  shown  in  Fig.  5.1  for  N  =  20,50 
and  =  0,  0.6.  (These  analytical  approximations  are  in  excellent  agree¬ 
ment  with  simulation  results.)  For  large  w,  divergence  of  the  expected 
estimator  from  the  actual  hazard  indicates  that  the  mean  square  error  in 
Eq.  5.11  is  due  mostly  to  bias;  hence,  poor  performance  of  the  present 
estimator  must  be  ascribed  to  the  inappropriate  choice  of  H,  which  contains 
no  model  close  to  the  actual  distribution,  especially  in  the  upper  tail. 

For  the  same  reason,  increasing  the  sample  size  above  N  =  20  produces  no 
appreciable  improvement. 

The  previous  analysis  is  readily  extended  to  the  case  of  nonexponential 
gaitma  families  (K  +  1).  The  maximum- likelihood  estimator  of  X  is  now 


X  =  K/W 


(5.12) 
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with  a  distribution  that  is  again  closely  approximated  by  an  inverted  gamma-1 
function.  Given  the  first  two  moments  of  W  or  the  parameters  Ky  and  Ay  that 
produce  the  best  gamma  approximation  to  Gy  (for  the  case  with  a£ne  =  0  these 
parameters  are  Ky  =  4.0  and  Ay  *  4.28),  given  the  shape  parameter  K  of  the 
K  family  and  the  sample  size  N,  one  can  use  tables  of  the  gamma  distribution 

/\  A  A 

to  calculate  fractiles  of  A  and  Gy(w).  For  example,  the  fractiles  of  Gy(w) 
in  Fig.  5.2  correspond  to  cr^  =  0,  K  =  3,  and  N  =  20.  Comparison  with 
Fig.  5.1  shows  considerable  reduction  of  bias.  Increasing  the  sample  size 
is  more  effective  than  for  the  exponential  family  because  a  larger  fraction 
of  the  mean  square  error  is  contributed  by  the  variance  term. 

b.  H  =  two-parameter  family  of  truncated-exponential  distributions 

The  likelihood  is  now  maximized  among  complementary  CDF  models  of  the 


Gy(w)  = 


1  -  Ml-e"15"'),  w  <  w. 


(5.13) 


w  >  W-, 


-ewi  -i 


in  which  8  and  w-|  are  unknown  positive  parameters  and  K-j  =  (1-e  )  .  The 

A 

maximum-likelihood  estimator  of  w-|  is  w.|  =  max  ,  but  no  simple  analytical 

A 

expression  is  available  for  B.  This  last  parameter  can  be  found  by  numeri¬ 
cally  solving  the  extreme  problem 


.max  /  e-eNW  } 

6  >  0  *•  ,  --Bwi  > 


(5.14) 


l.e‘^1 


Due  to  the  form  of  w-| ,  the  maximum-likelihood  estimator  of  JlnGy(w)  is  biased, 

A 

with  E[-Jln  Gy(w)]  =  »  for  all  w  and  all  finite  N.  It  is  therefore  more  mean¬ 
ingful  to  consider  sample  estimates  of  E[Gy(w)]  since  this  quantity  is  finite 
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nonzero  for  all  w.  Results  are  shown  as  curves  GTE  in  Figs.  5.3,  5.4,  and 
5.5.  Sample  estimates  of  the  standard  deviation  of  Gw(w)  and  sample  ranges 
are  also  indicated. 

In  maximum- likelihood  estimation,  the  one-parameter  exponential  and  the 
two-parameter  truncated  exponential  families  have  opposite  behavior  with 
regard  to  tail  probabilities.  The  exponential  family  (as  all  unbounded,  one- 
parameter  families)  extrapolates  the  upper  tail  from  the  body  of  the  distri¬ 
bution  where  the  data  lie,  whereas  the  truncated  exponential  family  (as  all 
families  with  upper- truncation  parameter)  makes  no  extrapolation  beyond  the 
largest  sample  value.  In  both  cases,  small-risk  estimates  are  strongly 
biased. 

c.  H  =  two-parameter  Wei  bull  family 

The  Weibull  distribution  with  parameters  a  and  c  has  a  complementary 
cumulative  function 

c 

Gw(w)  =  e"aw  ,w>0,  a>0,  c>0  (5.15) 

Given  the  independent  sample  W1 . WN,  the  maximum- likelihood  estimator 

of  c  can  be  calculated  numerically  from  the  condition 

N  c 

N  l  (W.)c  Jin  W. 

T  +  L  * n  W.  -  — - 7 - =  0  (5.16) 

c  1  =  1  V  \C 


and  then  a  can  be  found  from 


a  = 


I*  , 


(5.17) 


Probability  of  exc 


FIGURE  5.3  True  hazard  and  sample  statistics  of  maximum- 
likelihood  hazard  estimators  when  W  is  assumed 
to  have  truncated  exponential  (TE),  Wei  bull  (WE), 
or  lognormal  (LN)  distribution.  Twenty  samples  of 
size  20  from  the  true  distribution  of  W.  o.  »  0 
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FIGURE  5.4  True  hazard  and  sample  statistics  of  maximum- 
likelihood  hazard  estimators  when  W  is  assumed 
to  have  truncated  exponential  (TE),  Wei  bull  (WE), 
or  lognormal  (IN)  distribution.  Twenty  samples  of 

size  20  from  the  distribution  of  W.  o„  *  0.6. 
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Results  in  the  form  of  sample  mean  values  (curves  G^),  sample  standard 
deviations,  and  ranges  of  Gy(w)  are  shown  in  Figs.  5.3,  5.4,  and  5.5.  The 
Wei  bull  family  includes  distributions  that  are  close  to  the  time  distribu¬ 
tion  of  W.  Hence,  the  bias  is  smaller  and  the  percentage  reduction  of  mean 
square  error  from  increasing  the  sample  size  is  larger  than  for  distribution 
families  considered  before. 


d.  H  =  two-parameter  lognormal  family 

If  £n  W  is  assumed  to  have  normal  distribution  with  unknown  mean  m? 

2 

and  variance  ofcnW,  then  the  maximum-likelihood  estimators  of  these  parameters 
are 


xnW 


N 


tnW 


i  A, tn  “l 


i=l 

N 


(5.18) 


°tnW  =  .^i  Wi  '  m£nW^ 


and  the  maximum- likelihood  complementary  CDF  of  W  is 


£nW-mfnU 

Gw(w)  =  1  -  H—r — ) 


(5.19) 


3tn  W 


in  which  $  =  standard  normal  CDF.  Sample  results  are  shown  in  Figs.  5.3, 

a  i 

5.4,  and  5.5  as  G^  curves.  The  lognormal  estimator  performs  better  than 
the  exponential  but  not  as  well  as  the  Weibull  estimator,  mainly  because 

i 

the  tail  of  the  fitted  lognormal  distributions  falls  off  too  slowly.  j 


I 


103 


V 


5.3  Performance  of  Simple  Bayesian  Estimators 

Under  the  assumption  that  G^  belongs  to  a  parametric  family  K  of 
complementary  CDF's,  uncertainty  on  the  parameters  can  be  included  in 
seismic  hazard  calculations  through  Bayesian  analysis.  If  e  is  the 
vector  of  uncertain  parameters,  then  knowing  H  means  knowing  the  condi¬ 
tional  complementary  CDF,  GW|0(w,  _6).  The  prior  distribution  of  6_,  F0(jO» 
is  assumed  here  to  be  noninformative  of  either  of  the  following  types:  for 
components  with  unrestricted  values, 

dF*  ( 0) d 0  (5.20) 

ei 

whereas  for  components  that  are  bounded  on  one  side,  as  >  0q, 

dF^e>  - -5%  •  <>>e0  (5-21) 

A  priori  independence  is  assumed,  so  that  dF  (Q)  =  n  dF  (0.).  By  using  the 

—  i  °i 

distributions  in  Eqs.  5.20  and  5.21,  one  eliminates  dependence  of  estimated 
hazard  on  nonstatistical  information. 

Given  the  independent  sample  Z  =  {W-j ,  ...,  W^},  the  posterior  CDF  of  6. 
is,  from  Bayes1  theorem, 


dF”  (0)  a  dF^(0)  Lg|Z(0) 


in  which  the  likelihood  function  L0|Z  satisfies 

N 

L0|Z(_§)  a  .n  [-dGW|0(w,  ^)/dW  ] 

W-W^ 


(5.22) 


(5.23) 


# 
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Finally,  the  posterior  complementary  CDF  of  W  is  calculated  by  using  the 
Total  Probability  Theorem,  which  gives 


(w)  =  |  (w,  e)dF“  (©)  (5.24) 

This  last  function  is  the  Bayesian  estimator  of  G^(w).  More  explicit  results 
are  given  next  for  a  few  parametric  families  H. 

a.  K  =  one-parameter  exponential  or  two-parameter  truncated  exponential 
family 

Let  e  =  X  and  denote  by  GW|0  the  function  in  Eq.  5.6.  Then  Lx|Z(x) 

«  X  exp(-XNW)  with  W  =  sample  mean,  and  for  a  prior  distribution  of  X  in 
the  form  of  Eq.  5.21,  the  posterior  distribution  of  W  has  complementary 
CDF  (use  Eqs.  5.22  and  5.24), 


G'jj  (w)  -  (1  + 
w  NW 


-N 


,  w  >  0 


(5.25) 


It  can  be  shown  (Schumacher,  1977)  that  this  is  also  the  posterior  comple¬ 
mentary  CDF  of  W  if  H  is  the  family  of  truncated-exponential  distributions  in 

Eq.  5.13  and  X  and  w,  have  jointly  noninformative  prior  density,  f1  (x,w,) 

X,W1 

«  (xw-j)  ,  for  X  >  0,  w1  >  0. 


b.  H  =  two-parameter  Weibull  family 


For  e  =  [a,  c]  and  G^  in  Eq.  5. 1 5 ,  the  likelihood  function  of  given 


Z  is 


c,z(a,  c)«aNcN 


exp 


[  N  •,  N  „  , 

{-a  l  (W.)cl  n  (w.)c_1  (5.26) 

1  i=l  1  >  i=l  1 
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If  both  a  and  c  have  prior  distribution  in  the  form  of  Eq.  5.21,  then 
Eq.  5.24  yields 


«  r  •< 
n  (w,)c  1 

i=l  1 


r(N)dc 

~r  nr 

[  l  (w JC  +  wc] 
i=l  1 


(5.27) 


in  which  K  =  Gy(0)  is  a  normalization  constant.  The  integral  of  Eq.  5.27 
must  be  evaluated  numerically. 

c.  H  =  two-parameter  lognormal  family 
2 

Let  m£nW  and  have  prior  joint  density  function 


f;  2  (m,  o2)«X  ,  o2  *  0  (5.28) 

m£nW*  °i>nW  a 

Then  Bayesian  analysis  (Raiffa  and  Schlaifer,  1961)  shows  that  the  random 
variable 


a  n/h  i \  1/2  N  ao  1/2 

S(W)  =  (£n  W  -  m£nW)  [  Z  U"  wi  '  "W  ]  (5-29) 


has  posterior  Student-t  cumulative  distribution  Ft  ,  with  N-l  degrees  of 

CN-1 

freedom.  From  the  fact  that  S  is  a  monotonically  decreasing  function  of  W 
one  concludes, 


Gy(w)  =  P[S(W)  >  S(w)]  =  1  -  Ft  (S(w))  (5.30) 

w  XN-1 

Simulation  results  for  the  Bayesian  estimators  in  Eqs.  5.25,  5.27,  and 
5.30  are  shown  in  Figs.  5.6,  5.7,  and  5.8.  Previous  considerations  about 


FIGURE  5.6  True  hazard  and  sample  statistics  of  Bayesian  estimators 
when  W  is  assumed  to  have  exponential  or  truncated  ex¬ 
ponential  (EX),  Weibull  (WE),  or  lognormal  (LN)  distribu¬ 
tion.  Twenty  samples  of  size  20  from  the  true  distribution 


Probability  of  exceedance  given  an  event  with  a  *  10  cm/sec 


Probability  of  exceedance  given  an  event  with  a  >  10  cm/sec 
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FIGURE  5 


W  =  in  a  -  tnlO 


.8  True  hazard  and  sample  statistics  of  Bayesian  estimators 
when  W  Is  assumed  to  have  exponential  or  truncated 
exponential  (EX),  Welbull  (WE),  or  lognormal  (LN)  distribu¬ 
tion.  Seventeen  samples  of  size  50  from  the  true  distri¬ 
bution  of  W.  o.  =  0.6. 
line 


maximum-likelihood  estimators  still  apply.  However,  Bayesian  estimators 
with  noninformative  prior  distribution  of  the  parameters  tend  to  be  more 
conservative  in  the  range  of  high  ground  motion  intensities.  The  most 
evident  case  is  that  of  the  truncated-exponential  family,  but  differences 
in  mean  estimated  hazard  are  noticeable  also  in  the  Weibull  case  (compare 
Figs.  5.5  and  5.8).  One  reason  for  conservatism  is  that  Bayesian  estimation 
accounts  for  uncertainty  on  parameters  in  which  the  complementary  CDF  is 
very  nonlinear.  Bayesian  and  maximum-likelihood  estimators  coincide  in  the 
limit,  as  N  -*■  °°. 

5-4  More  Complicated  Estimators 

The  previous  results  indicate  that  hazard  estimators  are  very  sensitive 
to  the  choice  of  the  distribution  family,  Ji.  One  may  attempt  to  reduce  sen¬ 
sitivity  by  using  richer  families  of  distributions,  e.g.,  families  with 
three  or  four  parameters,  or  composite  distribution  sets  K  =  JK^,  ....  Hq| 
in  which  Ki  is  itself  a  distribution  family,  e.g.,  one  of  those  considered 
in  Secs.  4.2  and  4.3.  The  inference  rule  may  still  be  of  the  maximum- 
likelihood  or  Bayesian  type,  the  latter  in  the  form  of  a  weighted  average 
of  all  the  distributions  in  K.  A  Bayesian  approach  with  composite  H  has 
been  developed  by  Wood  et  al_.  (1974).  An  approximate,  simpler  "pseudo- 
Bayesian"  rule  has  been  suggested  by  Grigoriu  (1976),  of  the  type 

6W(w>  «  X  Pi’Li  GW.<w)  (5>31) 

in  which  P^'  =  prior  probability  that  G^  eK^. 

L.|  =  likelihood  of  the  maximum-likelihood  model  inK^. 

G"  =  posterior  Bayesian  complementary  CDF  of  W,  given  that 
wi 


no 


When  using  maximum- likelihood  or  pseudo-Bayesian  inference,  the  estimated 
hazard  function  critically  depends  on  the  likelihood  ratios  Lg/Lp  L3/L1 , 

....  Lq/I_i .  For  H  =  {K-j ,  K  g )  with  K-j  the  two-parameter  lognormal  family 
andl<2  the  two-parameter  Weibull  family,  and  for  10  independent  samples  of 
size  20  and  50,  ratios  L2/L^  are  given  in  Table  5.1.  All  samples  have  been 
pseudo-randomly  simulated  from  the  time  distribution  of  W  that  corresponds 
to  =  0.6.  Although  in  the  case  of  these  samples  a  Weibull  distribution 
always  exists  with  likelihood  higher  than  that  of  any  lognormal  distribution, 
the  large  sample  variability  of  Lg/l-i  implies,  at  least  for  N  =  20,  large 
statistical  variability  of  the  pseudo-Bayesian  estimator.  This  fact  in  turn 
indicates  that  recourse  to  a  richer  family  of  distributions  produces  better 
results  in  the  range  of  high  intensities,  only  if  the  sample  size  is  suf¬ 
ficiently  large  (e.g.,  N  >  50). 

All  model  estimators  M  (=  Gw)  and  model  families  JR  (-  M)  considered  up 
to  here  in  Sec.  5  are  of  the  site  type,  i.e.,  they  directly  describe  seismic 
hazard  at  the  site  of  interest.  If  historical  data  are  also  available  on 
earthquake  source  parameters  (e.g.,  epicentral  coordinates  and  magnitude), 
then  H  can  be  replaced  with  a  family  JR  of  regional  models,  e.g.,  of  the  same 
type  as  the  true  model  MT  in  Sec.  5.1.  As  an  example,  let  )R  collect 
all  regional  models  of  the  same  type  as  with  given  parameters  mg  =  4, 
b4  =  25,  a£ne  =  0  and  uncertain  parameters  mp  b,  r,  bp  b2,  and  bj.  The 
analyst  must  estimate  these  parameters  from  data  on  magnitude,  distance, 
and  peak  site  acceleration.  Statistical  sampling  is  from  the  true  mechanism 
with  parameters  in  Eq.  5.5  and  ofcng.  =0.6.  The  difference  between  assumed 
and  true  variance  of  the  attenuation  residual  term  e  is  therefore  the  only 
source  of  modeling  error.  Samples  have  the  form,  Z  =  { (M,  R,  a)p  .... 

(M,  R,  a)N>  in  which  (M,  R,  a ) ^  stands  for  magnitude,  epicentral  distance. 


Ill 


N 

■  20 

N  = 

50 

Sample 

# 

Vh 

■ 

L2/Ll 

Sample 

# 

l2/Li 

# 

L2/Ll 

1 

11.4 

6 

7.6 

1 

446 

6 

3354 

2 

4.1 

7 

13.7 

2 

4836 

7 

25 

3 

14.7 

8 

20.3 

3 

522 

8 

1615 

4 

1.3 

9 

403.2 

4 

2489 

44 

5 

354.1 

10 

9.8 

5 

22420 

H 

2161 

TABLE  5.1  Likelihood  ratio  Lp/L^  for  maxi mum- likelihood  Wei  bull  and 
lognormal  distributions.  Independent  samples  of  size 
N  are  simulated  from  the  true  distribution  of  W  for 
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and  peak  site  acceleration  for  the  ith  earthquake.  For  samples  of  this  type 
and  noninformative  prior  distribution 


fm1,b,r,£nb1,b2,b3  (™i.b,r,B,b2,b3)  “  m^5r 


(5.32) 


Schumacher  (1977)  found  an  expression  for  the  posterior  complementary  CDF, 
Gy(w).  The  average  of  this  function  over  20  independent  samples,  each  of 
size  N  =  20,  is  shown  in  Fig.  5.9  as  curve  G.  One-standard  deviation  in¬ 
tervals  are  also  indicated.  The  large  bias  is  due  primarily  to  uncertainty 
on  the  upper-bound  parameter  m-|  of  the  magnitude  distribution. 

The  most  important  consideration  from  the  error  analyses  of  this  section 
is  that,  when  hazard  estimates  are  based  solely  on  statistical  information, 
large  biases  may  result.  Biases  are  sometimes  due  to  the  estimation  pro¬ 
cedure  (recall  earlier  comparison  between  maximum-! ikeli hood  and  Bayesian 
inference  of  distributions  with  upper  truncation)  but  they  are  more  likely 
due  to  possible  errors  in  modeling  assumptions  (see  the  variation  of  mean 
hazard  with  K) . 

One  way  to  reduce  the  variability  of  hazard  estimates  is  to  utilize  non- 
statistical  information,  e.g.,  by  replacing  the  posterior  distribution  of  the 
parameters  in  Bayesian  analysis  with  subjectively  assessed  distributions. 

The  latter  distributions  are  usually  narrower  because  they  are  based  not 
only  on  historical  data,  but  also  on  nonstatistical  considerations  about 
local  geology,  earthquake  mechanisms,  and  similarity  with  other  seismic 
regions.  The  technique  relies  on  "expert  opinion"  and  is  becoming  in¬ 
creasingly  popular,  especially  for  the  assessment  of  hazard  in  regions  with 
poor  historical  records.  It  is  Itself  not  free  of  difficulties,  however,  as 
will  be  said  in  the  next  section. 


Probability  of  exceedance  given  an  event  with  a>  10  cm/sec' 


1 


2 


3 


4 


CSJ 


FIGURE  5.9  True  hazard  and  sample  statistics  (from  20  samples  of 
size  20)  of  a  Bayesian  estimator.  The  estimator 
corresponds  to  assuming  that  the  regional  earthquake 
mechanism  is  of  the  same  type  as  the  actual  mechanism, 
with  0£ne  =  0  and  a  few  uncertain  parameters.  For  the 
true  mechanism,  *  0.6. 
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6.  THE  USE  OF  EXPERT  OPINION 

The  basic  problem  when  turning  to  experts  for  the  quantification  of 
parameter  uncertainty  is  that  subjective  estimates  are  not  totally  reliable. 
For  this  reason  one  typically  asks  for  the  opinion  of  several  experts  and 
from  each  of  them  obtains  a  separate  hazard  estimate,  .  How  should  one 

A  A 

interpret  the  difference  between  H.  and  H.?  What  is  the  procedure  to  com- 
bine  the  H^  and  produce  a  single  (sometimes  called  consensus)  hazard  esti¬ 
mate  $?  How  does  accuracy  of  H  depend  on  the  number  n  of  experts?  On 
their  degree  of  expertise? 

The  first  step  toward  answering  these  questions  is  to  make  a  concep¬ 
tual  model  of  professional  judgement.  Conceptual  models  should  then  be 
given  mathematical  form  and  analyzed.  The  present  section  follows  this 
ideal  path  by  first  introducing  two  alternative  interpretations  of  the 
cause  for  diversity  in  expert  opinion.  Each  interpretation  contains  ele¬ 
ments  of  truth.  The  mathematical  model  that  follows  in  Sec.  6.2  incor¬ 
porates  features  of  both  interpretations.  It  is  flexible  and  simple. 

There  is  obviously  no  contention  of  completeness  in  dealing  with  a  topic 
which  has  been  under  study  for  decades,  primarily  by  nonengineers;  see 
for  example  Winkler  (1968),  de  Groot  (1974),  Morris  (1974),  Dal  key  (1975). 

6.1  Conceptual  Models 

There  are  two  seemingly  conflicting  interpretations  of  the  diversity 
of  expert  opinion.  One  interpretation  views  each  estimate  as  a  noisy 
observation  of  the  true  hazard  Hj-,  the  other  considers  each  estimate  as  a 
sample  from  a  random  process  that  expresses  objective  professional 
uncertainty  about  Hy. 
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Observation  model.  Through  mathematical  and  cognitive  procedures, 
each  expert  operates  on  statistical  and  physical  information  to  produce  an 
estimate  of  Hy.  The  difference  between  and  Hy  is  attributed  to 
statistical  estimation  or  observation  errors,  as  if  experts  were  "ob¬ 
serving"  the  true  hazard  by  means  of  imperfect  instruments.  Different 
experts  may  use  different  sets  of  data  and  different  theories  and  pro¬ 
cedures,  thus  producing  different  estimates. 

Sampling  model.  It  is  recognized  at  the  outset  that  statistical  infor¬ 
mation  and  physical  knowledge  are  limited  so  that  no  matter  how  many  experts 

are  questioned,  there  will  always  be  uncertainty  about  the  actual  hazard. 

/\ 

Estimates  Hi  are  considered  as  random  samples  that  reflect  this  irreducible, 
professional  uncertainty.  Therefore,  one  should  not  reduce  the  to  a 
single  "consensus"  estimate,  H.  Rather,  the  full  diversity  of  opinions 
should  be  retained  and  used  to  quantify  professional  uncertainty  on  Hy. 

Both  models  are  valuable.  For  example,  it  is  logical  to  think  of  an 
objective  irreducible  uncertainty  on  Hy  as  predicated  by  the  sampling  model. 
However,  viewing  expert  estimates  as  I.I.D.  variables  from  a  hypothetical 
"professional  distribution"  is  not  correct.  Not  only  does  each  expert  oper¬ 
ate  on  a  different  fraction  of  the  total  information  available,  but  he  also 
makes  errors  depending  on  professional  experience,  degree  of  familiarity 
with  hazard  models,  etc.  One  might  think,  therefore,  of  some  additional 
type  of  error  in  the  H^,  of  the  type  emphasized  by  the  observational  model. 
This  will  be  called  personal  error.  Whereas  professional  uncertainty 
is  the  same  for  all  experts,  personal  uncertainty  varies.  This  varia¬ 
tion  is  one  of  the  reasons  why  one  may  want  to  give  different  weights 
to  different  estimates.  In  the  next  section,  mathematical  models  will 
be  developed  to  represent  all  these  components  of  uncertainty.  First  a 
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quantitative  model  of  the  observational  point  of  view  will  be  developed 
and  then  it  will  be  generalized  to  incorporate  the  idea  of  professional 
uncertainty.  It  is  always  assumed  that  the  objective  is  to  estimate 
(uncertainty  on)  the  mean  rate  A  at  which  a  given  ground  acceleration 
level  a  is  exceeded  at  the  site.  The  estimate  of  A  based  on  information 
from  the  ith  expert  is  denoted  by  a.. . 

6.2  Mathematical  Models 

According  to  the  observation  model,  the  ith  log  hazard  estimate,  «.nA^ 
is  considered  to  be  a  noisy  observation  of  fin  A  with  random  error  ei  = 
fin  A.j  -  fin  A,  so  that 


fin  A^  =  fin  A  +  (6.1) 

The  fact  that  experts  use  basically  the  same  historical  data,  geological 
information,  seismicity  theories,  and  hazard  procedures  makes  the  error 
terms  in  Eq.  6.1  probabilistically  dependent.  Suppose  for  convenience 
that  the  vector  £  =  ....  en]T  has  multivariate  normal  distribution 

with  zero  mean  and  covariance  matrix  .  If,  before  obtaining  any  expert 
estimate,  the  hazard  fin  A  is  independent  of  £  and  has  noninformative  flat 
prior  distribution,  then  a  posteriori  fin  A  has  normal  distribution  with  mean 

m£nA  =  (1T  in  A)/(1T  1)  (6.2) 


and  variance 


-  V(1T 


1) 


(6.3) 
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where  1  denotes  the  column  vector  with  n  unit  components  and  Hn  A  is  the 

a  a  /n  t 

vector  of  estimates,  Jin  A  =  [Jin  . in  An]  .  In  the  case  of  only  one 

expert,  this  posterior  distribution  becomes  simply  N(jln  A,  o2). 

Due  to  linearity  of  the  observation  model  in  Eq.  6.1  and  normality  of 
£,  the  posterior  mean  value  in  Eq.  6.2  is  a  linear  combination  of  the  ex¬ 
pert  estimates,  with  coefficients  that  depend  on  the  covariance  matrix  of 
the  errors.  The  posterior  variance,  Eq.  6.3,  is  a  nonincreasing  function 
of  the  number  of  experts  and  for  given  n,  increases  with  the  variances  and 
the  correlation  coefficients  of  the  errors. 

If  the  errors  were  independent  and  identically  distributed  so  that 
2 

=  o£  _I,  then  one  would  find  that  a  posteriori 


Jin  A  ~  N  (1  T  Jtn  A.  ,  1  o2) 

n  >i  i  n  e 


(6.4) 


2 

The  posterior  variance  of  in  A  is  larger  than  o£/n  if  errors  are  positively 
correlated.  A  case  with  correlation  which  leads  to  explicit  results  is  that 
in  which  experts  are  treated  as  equally  credible  and  correlation  between 
estimation  errors  is  the  same,  irrespective  of  the  identity  of  the  experts. 
Covariance  matrices  that  express  these  conditions  have  the  equicorrelated 
form. 


l  = 
— z 


e 


(6.5) 


This  last  model  may  not  be  accurate  in  practice  if  the  identity  of  the 
experts  is  known,  but  is  realistic  in  the  case  of  randomization,  i.e.. 
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when  pairing  between  experts  and  estimates  is  not  known  or  is  ignored. 

The  technique  of  randomization  produces  suboptimal  estimates  because  it 
neglects  available  information,  but  has  the  benefit  of  reducing  the  com¬ 
plexity  of  analysis  and  the  number  of  parameters  (to  only  n,  a2,  and  p).  For 
in  Eq.  6.5,  the  posterior  distribution  of  An  A  is 

n 

in  A~N  (1  J  An  ^  ,  a\  1-LlQ-D.P.)  (6.6) 

The  quantity  ]  +  )p  »  which  can  be  interpreted  as  the  equivalent  number 
of  independent  experts,  is  plotted  in  Fig.  6.1  as  a  function  of  n  for  selected 
values  of  the  correlation  coefficient  p.  Notice  that,  as  n  +  the  posterior 
variance  tends  asymptotically  to 

°l  =  l™  °*"A  =  P  ^  (6‘7) 

which  is  a  quantity  larger  than  zero,  unless  p  =  0.  This  means  that  an 

infinite  number  of  expert  estimates  with  common  correlation  coefficient 

p  has  the  same  information  content  as  a  single  estimate  with  error  vari- 
2 

ance  pa£. 

A  possible  interpretation  of  the  asymptotic  posterior  distribution 
2 

N(moo,  oj  is  that  this  distribution  represents  irreducible  professional 
uncertainty  about  the  true  seismic  hazard.  This  is  the  uncertainty  after 
considering  the  opinion  of  all  living  experts;  it  cannot  be  reduced  because 
it  is  due  to  objective  limitation  of  information  and  knowledge. 

Let  Cq  be  a  random  variable  with  the  asymptotic  distribution  of  An  A 
and  define  variables  n-j ,  ...,  with  distribution  r^.  -  N(0,  o  ).  The 


1 


2 
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FIGURE  6. 


3  4  5  6 

Number  of  experts,  n 


1  Equivalent  number  of  independent 
experts  as  a  function  of  n  and  p 
(equicorrelated  model). 
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quantity  ni  can  be  regarded  as  the  component  of  personal  estimation  error 
for  the  ith  expert.  With  Cq  and  the  mutually  independent  variables, 
the  model  of  Eq.  6.1  can  be  written 

£n  A^  =  £n  A  +  e_  +  (6.8) 

in  which  eg  +  ni  has  been  substituted  for  c^.  For  example,  the  covariance 
matrix  in  Eq.  6.5  results  from  setting  =  (l-p)o£  for  all  experts. 

It  is  clear  from  Eq.  6.8  that  the  present  mathematical  model  considers 
professional  uncertainty  as  systematic  error  (bias)  and  associates  it  with 
the  random  variable  Cg.  On  the  contrary,  the  sampling  viewpoint  predicates 
that  professional  uncertainty  is  responsible  for  variation  among  expert 
estimates,  not  for  bias.  It  is  believed  that  both  concepts  are  true  and 
that  professional  uncertainty  has  in  fact  both  a  bias  and  a  random  component-- 
bias  due  to  limitation  of  information,  and  random  component  caused  by  alter¬ 
native  possible  interpretations  of  seismicity  and  alternative  models  for 
hazard  calculation.  The  model  of  Eq.  6.8  should  then  be  generalized  to 

In  =  £n  A  +  £q  +  £.  +  r^.  (6.9) 


in  which 

2 

£q  -  N(0,  Oq)  is  random  bias  from  limited  information. 

p 

£.  _  N(0,  oJ  is  a  random  component  associated  with  uncertainty 
1  ^  on  the  interpretation  of  seismicity  and  on  the 

model  for  hazard  calculation. 

2 

ni  -  N(0,  o^)  is  personal  error. 
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Correlation  may  exist  among  some  of  these  variables.  Professional  uncer- 

2  2 

tainty  corresponds  to  the  random  variable,  y  =  £g  +  ei  -  N(0,  Og  +  o^). 

The  objective  of  analysis  is  to  estimate  the  distribution  of  y,  j,  e.  ,  the 
2  2 

variance  terms  Og  and  a^. 

Unfortunately,  it  is  clear  that  the  model  postulated  by  Eq.  6.9  is  not 

A  ss 

inferable  from  the  assumed  data  structure,  (£n  A-| ,  ....  Jin  An>.  Specifically, 
the  variance  of  eg  cannot  be  estimated  and  it  is  ii, .possible  to  separate 
n-  from  r.j  unless  correlation  is  introduced.  One  can  think  of  other, 
realistic  data  structures  that  would  make  estimation  feasible,  e.g.,  struc¬ 
tures  that  are  typical  of  analysis-of-variance  problems. 

Let  us  start  with  separation  between  £;.  and  nn- .  A  way  to  estimate  the 
2  2 

variances  o ^  and  is  to  first  specify  a  realistic  set  of  alternative 
seismicity  interpretations  and  hazard  models  and  form  a  matrix  of  estimates 

/\ 

Hn  A.  =  estimate  of  fcn  A  using  information 
j  from  expert  i  and  the  jth  combination 
of  seismicity  interpretation  and 
hazard  model 

(6.10) 

Standard  methods  of  analysis  of  variance  (e.g..  Draper  and  Smith,  1966)  will 

2  2 

then  provide  estimates  of  or  and  o  . 

K  n 

2 

The  estimation  of  Og  is  somewhat  more  complicated.  One  way  to  obtain  in- 

2 

formation  on  Og  is  through  analyses  of  the  type  in  Sec.  5,  i.e.,  by  pseudo- 
randomly  generating  alternative  data  sets  (possibly  also  alternative  seismic 
environments)  and  by  calculating  the  associated  variation  of  estimated  hazard. 
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6.3  Actual  Variability  of  Expert  Opinion 

In  a  recent  study  of  Seismic  Hazard  in  the  Eastern  United  States 
(TERA,  1979),  10  seismologists  responded  to  a  comprehensive  questionnaire 
which  included  information  on  source  geometry,  seismicity  rate  and  magni¬ 
tude  distribution  for  each  source,  and  form  and  parameters  of  the  attenu¬ 
ation  function.  Before  answering  the  questionnaire,  experts  were  provided 
with  available  maps  of  seismicity  and  with  exhaustive  catalogues  of  his¬ 
torical  events.  In  spite  of  this  attempt  to  equalize  the  information 
available  to  all  experts,  opinions  ranged  widely  on  a  few  critical  param¬ 
eters,  notably,  on  the  appropriate  configuration  of  earthquake  sources,  on 
the  slope  of  the  recurrence  relationship,  and  on  the  upper  bound  size. 

Based  on  response  from  each  expert,  10  separate  estimates  of  hazard  at  vari 
ous  locations  of  the  eastern  and  central  United  States  were  obtained.  Re¬ 
sults  at  one  site,  in  terms  of  pseudo-velocity  spectra  for  a  1000-year  re¬ 
turn  period,  are  shown  in  Fig.  6.2.  The  range  of  estimates  corresponds  to 
factors  on  the  spectrum  of  typically  2  or  3,  with  peak  values  of  about  7 
for  long-period  oscillators.  With  reference  to  the  model  of  Eq.  6.9,  this 
variability  should  be  attributed  mainly  to  the  term  and  less  to  the 
term  ^  because  of  the  similarity  of  hazard  analysis  procedures  (but  the 
interpretation  of  seismicity  was  allowed  to  vary  from  expert  to  expert). 

A  self-ranking  method  (experts  were  asked  to  quantify  their  degree 
of  confidence  in  answering  various  parts  of  the  questionnaire)  was  devised 

/V 

to  finally  produce  a  single  hazard  estimate  H  as  a  weighted  average  of  the 

A 

individual  curves.  However,  precision  of  H  was  not  quantified. 

Another  example  is  shown  in  Fig.  6.3  (from  Cornell,  1980).  In  this 


case,  three  independent  consulting  teams  obtained  estimates  of  seismic 
hazard  at  the  site  of  the  Diablo  Canyon  nuclear  power  plant,  in  California. 


Ril  EXPERTS  —  CONNECTICUT  TRNKEE  —  iOOO  TERR  RETURN  PERIOD 
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FIGURE  6.2  Diversity  of  expert  opinion  in  the 
assessment  of  seismic  hazard  at 
Connecticut  Yankee  power  station 
site  (from  TERA,  1979) 
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Although  the  sample  is  very  small,  the  range  of  estimates  corresponds  to 
a  factor  of  about  10  on  probability  given  a  and  of  about  3  on  acceleration 
given  the  probability  level. 

Another  useful  source  of  information  on  the  bias  and  error  of  sub¬ 
jective  quantification  and  updating  of  uncertainty  is  the  data  collected 
by  experimental  psychologists  (Phillips,  Mays,  and  Edwards,  1966;  Winkler, 
1967;  Vlek,  1973;  Stael  Von  Holstein,  1971;  Rapoport  and  Wallsten,  1972; 
to  name  but  a  few).  Distortions  due  to  conservatism  and  hedging  are  rather 
well  documented;  so  are  techniques  for  controlling  them,  such  as  feedback 
and  cross-checking. 


7.  CONCLUSIONS 


It  is  just  natural  that  first  theories  and  models  are  developed  and 
then  their  limitations  are  understood.  In  the  case  of  seismic  hazard 
analysis,  the  first  step  is  essentially  completed  but  not  the  second, 
so  that  quoting  from  a  contemporary  survey  of  the  field  (Cornell,  1980), 
"the  probabilistic  procedures  of  applied  seismic  hazard  analysis  are  at 
the  moment  very  stable." 

By  contrast,  understanding  of  the  accuracy  of  even  the  standard 
methods  is  at  best  incomplete.  The  very  definition  of  accuracy  (with 
respect  to  what?)  has  never  been  given  in  clear,  unequivocal  terms. 
Studies  of  the  subject  are  fragmentary  and  uncoordinated;  they  look 
at  different  sources  of  error  and  use  different  methods  and  language 
in  a  rather  unscientific  fashion.  Besides  the  conceptual  need  to  give 
order  to  the  subject,  there  is  the  practical  necessity  of  deciding  what 
to  do  about  uncertainty  on  the  true  hazard.  While  waiting  for  an 
answer,  there  is  a  consensus,  de  facto,  that  one  should  be  conservative, 
although  nobody  has  come  up  with  a  measure  of  conservatism,  even  less 
with  a  rationale. 

The  purpose  of  this  report  is  to  gather  available  information  on 
errors  in  seismic  hazard  analysis  and,  more  important,  to  provide  a 
conceptual  framework  to  the  subject.  Toward  the  latter  end,  it  was 
found  necessary  to  first  define  a  "true  probabilistic  hazard"  against 
which  accuracy  could  be  measured.  An  objective  definition  of  hazard  has 
been  given  In  Sec.  2,  as  the  relative  frequency  in  time  with  which  earth¬ 
quake  events  of  engineering  Interest  occur  at  the  site.  The  body  of 
the  report— Secs.  4,  5,  and  6— deals  with  errors  in  the  estimation 
of  true  hazard— statistical  errors  In  Secs.  4  and  5,  judgmental 


errors  in  Sec.  6.  Specifically,  Sec.  4  looks  at  errors  from  modeled 
uncertainty  (on  parameters  or  on  alternative  hazard  procedures)  whereas 
Sec.  5  deals  with  unmodeled  errors  from  making  wrong  assumptions  about 
the  relevant  laws  of  nature  (form  of  the  true  hazard  function  or  charac¬ 
teristics  of  seismicity).  Errors  of  the  latter  type  are  more  elusive  and 
of ten  more  severe  than  those  that  are  explicitly  modeled.  Hence,  judging 
accuracy  only  from  modeled  uncertainty  (e.g.,  about  parameters  such  as 
the  Richter  b  slope  or  the  upper  bound  magnitude)  may  give  a  false  sense 
of  confidence  in  a  model. 

With  the  increasing  use  of  nonstatistical  information  in  the  form  of 
professional  judgment,  the  study  of  errors  in  the  subjective  estimation 
and  updating  of  uncertainty  becomes  part  of  the  accuracy  problem.  A 
related  issue  is  that  of  how  to  combine  hazard  estimates  from  different 
seismologist  experts.  Reflections  on  this  problem  are  offered  in  Sec.  6. 

It  is  realized  that  the  engineering  accuracy  and  value  of  a  model  can¬ 
not  be  quantified  only  in  terms  of  errors.  Consequences  of  errors  rather 
than  errors  are  important,  uncertainty  on  seismic  risk  rather  than  on 
seismic  hazard.  Although  most  of  the  literature  to  date  deals  with  the 
latter  topic,  a  section  (Sec.  3)  has  been  included  here  on  seismic  risk 
and  its  dependence  on  parameter  uncertainty. 

It  is  apparent  that  more  research  is  needed  on  virtually  every  aspect 
of  this  complex  problem  and  that  it  is  of  utmost  importance  that  it  be 
carried  out  with  much  more  coordination  than  we  have  seen  in  the  past. 
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